package MathCore;

/**
 * <b>Java integer implementation of 63-bit precision floating point.</b>
 * <br><i>Version 1.09</i>
 *
 * <p>Copyright 2003-2006 Roar Lauritzsen <roarl@pvv.org>
 *
 * <blockquote>
 * <p>This library is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * <p>This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * <p>The following link provides a copy of the GNU General Public License:
 * <br>&nbsp;&nbsp;&nbsp;&nbsp;<a
 * href="http://www.gnu.org/licenses/gpl.txt">http://www.gnu.org/licenses/gpl.txt</a>
 * <br>If you are unable to obtain the copy from this address, write to
 * the Free Software Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA 02111-1307 USA
 * </blockquote>
 *
 * <p><b>General notes</b>
 * <ul>
 *
 *   <li><code>Real</code> objects are not immutable, like Java
 *   <code>Double</code> or <code>BigDecimal</code>. This means that
 *   you should not think of a <code>Real</code> object as a "number",
 *   but rather as a "register holding a number". This design choice
 *   is done to encourage object reuse and limit garbage production
 *   for more efficient execution on e.g. a limited MIDP device. The
 *   design choice is reflected in the API, where an operation like
 *   {@link #add(Real) add} does not return a new object containing
 *   the result (as with {@link
 *   java.math.BigDecimal#add(java.math.BigDecimal) BigDecimal}), but
 *   rather adds the argument to the object itself, and returns
 *   nothing.
 *
 *   <li>This library implements infinities and NaN (Not-a-Number)
 *   following the IEEE 754 logic. If an operation produces a result
 *   larger (in magnitude) than the largest representable number, a
 *   value representing positive or negative infinity is generated. If
 *   an operation produces a result smaller than the smallest
 *   representable number, a positive or negative zero is
 *   generated. If an operation is undefined, a NaN value is
 *   produced. Abnormal numbers are often fine to use in further
 *   calculations. In most cases where the final result would be
 *   meaningful, abnormal numbers accomplish this,
 *   e.g. atan(1/0)=&pi;/2. In most cases where the final result is
 *   not meaningful, a NaN will be produced.
 *
 *   <li>Error bounds listed under <a href="#method_detail">Method
 *   Detail</a> are calculated using William Rossi's <a
 *   href="http://dfp.sourceforge.net/">rossi.dfp.dfp</a> at 40
 *   decimal digits accuracy. Error bounds may increase when results
 *   approach zero or infinity. The abbreviation {@link
 *   Math#ulp(double) ULP} means Unit in the Last Place. An error
 *   bound of � ULP means that the result is correctly rounded. The
 *   relative execution time listed under each method is the average
 *   from running on SonyEricsson T610 (R3C), K700i, and Nokia 6230i.
 *
 *   <li>The library is not thread-safe. Static <code>Real</code>
 *   objects are used extensively as temporary values to avoid garbage
 *   production and the overhead of <code>new</code>. To make the
 *   library thread-safe, references to all these static objects must
 *   be replaced with code that instead allocates new
 *   <code>Real</code> objects in their place.
 *
 *   <li>There is one bug that occurs again and again and is really
 *   difficult to debug. Although the pre-calculated constants are
 *   declared <code>static final</code>, Java cannot really protect
 *   the contents of the objects in the same way as
 *   <code>const</code>s are protected in C/C++. Consequently, you can
 *   accidentally change these values if you send them into a function
 *   that modifies its arguments. If you were to modify {@link #ONE
 *   Real.ONE} for instance, many of the succeeding calculations would
 *   be wrong because the same variable is used extensively in the
 *   internal calculations of Real.java.
 *
 * </ul>
 */
public final class Real
{
  /**
   * The mantissa of a <code>Real</code>. <i>To maintain numbers in a
   * normalized state and to preserve the integrity of abnormal
   * numbers, it is discouraged to modify the inner representation of
   * a <code>Real</code> directly.</i>
   *
   * <p>The number represented by a <code>Real</code> equals:<br>
   * &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;-1<sup>sign</sup>&nbsp;�&nbsp;mantissa&nbsp;�&nbsp;2<sup>-62</sup>&nbsp;�&nbsp;2<sup>exponent-0x40000000</sup>
   *
   * <p>The normalized mantissa of a finite <code>Real</code> must be
   * between <code>0x4000000000000000L</code> and
   * <code>0x7fffffffffffffffL</code>. Using a denormalized
   * <code>Real</code> in <u>any</u> operation other than {@link
   * #normalize()} may produce undefined results. The mantissa of zero
   * and of an infinite value is <code>0x0000000000000000L</code>.
   *
   * <p>The mantissa of a NaN is any nonzero value. However, it is
   * recommended to use the value
   * <code>0x4000000000000000L</code>. Any other values are reserved
   * for future extensions.
   */
  public long mantissa;

  /**
   * The exponent of a <code>Real</code>. <i>To maintain numbers in a
   * normalized state and to preserve the integrity of abnormal
   * numbers, it is discouraged to modify the inner representation of
   * a <code>Real</code> directly.</i>
   *
   * <p>The exponent of a finite <code>Real</code> must be between
   * <code>0x00000000</code> and <code>0x7fffffff</code>. The exponent
   * of zero <code>0x00000000</code>.
   *
   * <p>The exponent of an infinite value and of a NaN is any negative
   * value. However, it is recommended to use the value
   * <code>0x80000000</code>. Any other values are reserved for future
   * extensions.
   */
  public int exponent;

  /**
   * The sign of a <code>Real</code>. <i>To maintain numbers in a
   * normalized state and to preserve the integrity of abnormal
   * numbers, it is discouraged to modify the inner representation of
   * a <code>Real</code> directly.</i>
   *
   * <p>The sign of a finite, zero or infinite <code>Real</code> is 0
   * for positive values and 1 for negative values. Any other values
   * may produce undefined results.
   *
   * <p>The sign of a NaN is ignored. However, it is recommended to
   * use the value <code>0</code>. Any other values are reserved for
   * future extensions.
   */
  public byte sign;

  /**
   * Set to <code>false</code> during numerical algorithms to favor
   * accuracy over prettyness. This flag is initially set to
   * <code>true</code>.
   *
   * <p>The flag controls the operation of a subtraction of two
   * almost-identical numbers that differ only in the last three bits
   * of the mantissa. With this flag enabled, the result of such a
   * subtraction is rounded down to zero. Probabilistically, this is
   * the correct course of action in an overwhelmingly large
   * percentage of calculations.
   *
   * <p>However, certain numerical algorithms such as differentiation
   * depend on keeping maximum accuracy during subtraction.
   */
  public static boolean magicRounding = true;

  /**
   * A <code>Real</code> constant holding the exact value of 0.
   * Among other uses, this value is used as a result when a positive
   * underflow occurs.
   */
  public static final Real ZERO = new Real(0,0x00000000,0x0000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 1.
   */
  public static final Real ONE = new Real(0,0x40000000,0x4000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 2.
   */
  public static final Real TWO = new Real(0,0x40000001,0x4000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 3.
   */
  public static final Real THREE = new Real(0,0x40000001,0x6000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 5.
   */
  public static final Real FIVE = new Real(0,0x40000002,0x5000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 10.
   */
  public static final Real TEN = new Real(0,0x40000003,0x5000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 100.
   */
  public static final Real HUNDRED= new Real(0,0x40000006,0x6400000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of 1/2.
   */
  public static final Real HALF = new Real(0,0x3fffffff,0x4000000000000000L);
  /**
   * A <code>Real</code> constant that is closer than any other to 1/3.
   */
  public static final Real THIRD = new Real(0,0x3ffffffe,0x5555555555555555L);
  /**
   * A <code>Real</code> constant that is closer than any other to 1/10.
   */
  public static final Real TENTH = new Real(0,0x3ffffffc,0x6666666666666666L);
  /**
   * A <code>Real</code> constant that is closer than any other to 1/100.
   */
  public static final Real PERCENT= new Real(0,0x3ffffff9,0x51eb851eb851eb85L);
  /**
   * A <code>Real</code> constant that is closer than any other to the square
   * root of 2.
   */
  public static final Real SQRT2 = new Real(0,0x40000000,0x5a827999fcef3242L);
  /**
   * A <code>Real</code> constant that is closer than any other to the square
   * root of 1/2.
   */
  public static final Real SQRT1_2= new Real(0,0x3fffffff,0x5a827999fcef3242L);
  /**
   * A <code>Real</code> constant that is closer than any other to 2&pi;.
   */
  public static final Real PI2 = new Real(0,0x40000002,0x6487ed5110b4611aL);
  /**
   * A <code>Real</code> constant that is closer than any other to &pi;,
   * the ratio of the circumference of a circle to its diameter.
   */
  public static final Real PI = new Real(0,0x40000001,0x6487ed5110b4611aL);
  /**
   * A <code>Real</code> constant that is closer than any other to &pi;/2.
   */
  public static final Real PI_2 = new Real(0,0x40000000,0x6487ed5110b4611aL);
  /**
   * A <code>Real</code> constant that is closer than any other to &pi;/4.
   */
  public static final Real PI_4 = new Real(0,0x3fffffff,0x6487ed5110b4611aL);
  /**
   * A <code>Real</code> constant that is closer than any other to &pi;/8.
   */
  public static final Real PI_8 = new Real(0,0x3ffffffe,0x6487ed5110b4611aL);
  /**
   * A <code>Real</code> constant that is closer than any other to <i>e</i>,
   * the base of the natural logarithms.
   */
  public static final Real E = new Real(0,0x40000001,0x56fc2a2c515da54dL);
  /**
   * A <code>Real</code> constant that is closer than any other to the
   * natural logarithm of 2.
   */
  public static final Real LN2 = new Real(0,0x3fffffff,0x58b90bfbe8e7bcd6L);
  /**
   * A <code>Real</code> constant that is closer than any other to the
   * natural logarithm of 10.
   */
  public static final Real LN10 = new Real(0,0x40000001,0x49aec6eed554560bL);
  /**
   * A <code>Real</code> constant that is closer than any other to the
   * base-2 logarithm of <i>e</i>.
   */
  public static final Real LOG2E = new Real(0,0x40000000,0x5c551d94ae0bf85eL);
  /**
   * A <code>Real</code> constant that is closer than any other to the
   * base-10 logarithm of <i>e</i>.
   */
  public static final Real LOG10E = new Real(0,0x3ffffffe,0x6f2dec549b9438cbL);
  /**
   * A <code>Real</code> constant holding the maximum non-infinite positive
   * number = 4.197e323228496.
   */
  public static final Real MAX = new Real(0,0x7fffffff,0x7fffffffffffffffL);
  /**
   * A <code>Real</code> constant holding the minimum non-zero positive
   * number = 2.383e-323228497.
   */
  public static final Real MIN = new Real(0,0x00000000,0x4000000000000000L);
  /**
   * A <code>Real</code> constant holding the value of NaN (not-a-number).
   * This value is always used as a result to signal an invalid operation.
   */
  public static final Real NAN = new Real(0,0x80000000,0x4000000000000000L);
  /**
   * A <code>Real</code> constant holding the value of positive infinity.
   * This value is always used as a result to signal a positive overflow.
   */
  public static final Real INF = new Real(0,0x80000000,0x0000000000000000L);
  /**
   * A <code>Real</code> constant holding the value of negative infinity.
   * This value is always used as a result to signal a negative overflow.
   */
  public static final Real INF_N = new Real(1,0x80000000,0x0000000000000000L);
  /**
   * A <code>Real</code> constant holding the value of negative zero.
   * This value is used as a result e.g. when a negative underflow occurs.
   */
  public static final Real ZERO_N = new Real(1,0x00000000,0x0000000000000000L);
  /**
   * A <code>Real</code> constant holding the exact value of -1.
   */
  public static final Real ONE_N = new Real(1,0x40000000,0x4000000000000000L);

  private static final int clz_magic = 0x7c4acdd;
  private static final byte[] clz_tab =
    { 31,22,30,21,18,10,29, 2,20,17,15,13, 9, 6,28, 1,
      23,19,11, 3,16,14, 7,24,12, 4, 8,25, 5,26,27, 0 };

  /**
   * Creates a new <code>Real</code> with a value of zero.
   */
  public Real() {
  }

  /**
   * Creates a new <code>Real</code>, assigning the value of another
   * <code>Real</code>. See {@link #assign(Real)}.
   *
   * @param a the <code>Real</code> to assign.
   */
  public Real(Real a) {
    { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
  }

  /**
   * Creates a new <code>Real</code>, assigning the value of an integer.
   * See {@link #assign(int)}.
   *
   * @param a the <code>int</code> to assign.
   */
  public Real(int a) {
    assign(a);
  }

  /**
   * Creates a new <code>Real</code>, assigning the value of a long integer.
   * See {@link #assign(long)}.
   *
   * @param a the <code>long</code> to assign.
   */
  public Real(long a) {
    assign(a);
  }

  /**
   * Creates a new <code>Real</code>, assigning the value encoded in a
   * <code>String</code> using base-10.
   * See {@link #assign(String)}.
   *
   * @param a the <code>String</code> to assign.
   */
  public Real(String a) {
    assign(a,10);
  }

  /**
   * Creates a new <code>Real</code>, assigning the value encoded in a
   * <code>String</code> using the specified number base.
   * See {@link #assign(String,int)}.
   *
   * @param a the <code>String</code> to assign.
   * @param base the number base of <code>a</code>. Valid base values are
   *     2, 8, 10 and 16.
   */
  public Real(String a, int base) {
    assign(a,base);
  }

  /**
   * Creates a new <code>Real</code>, assigning a value by directly setting the
   * fields of the internal representation. The arguments must
   * represent a valid, normalized <code>Real</code>. This is the
   * fastest way of creating a constant value.
   * See {@link #assign(int,int,long)}.
   *
   * @param s {@link #sign} bit, 0 for positive sign, 1 for negative sign
   * @param e {@link #exponent}
   * @param m {@link #mantissa}
   */
  public Real(int s, int e, long m) {
    { this.sign=(byte)s; this.exponent=e; this.mantissa=m; };
  }

  /**
   * Creates a new <code>Real</code>, assigning the value previously encoded
   * into twelve consecutive bytes in a byte array using
   * {@link #toBytes(byte[],int) toBytes}.
   * See {@link #assign(byte[],int)}.
   *
   * @param data byte array to decode into this <code>Real</code>.
   * @param offset offset to start encoding from. The bytes
   *     <code>data[offset]...data[offset+11]</code> will be
   *     read.
   */
  public Real(byte [] data, int offset) {
    assign(data,offset);
  }

  /**
   * Assigns this <code>Real</code> the value of another <code>Real</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to assign.
   */
  public void assign(Real a) {
    if (a == null) {
      makeZero();
      return;
    }
    sign = a.sign;
    exponent = a.exponent;
    mantissa = a.mantissa;
  }

  /**
   * Assigns this <code>Real</code> the value of an integer.
   * All integer values can be represented without loss of accuracy.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = (double)a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.6
   * </td></tr></table>
   *
   * @param a the <code>int</code> to assign.
   */
  public Real assign(int a) {
    if (a==0) {
      makeZero();
      return this;
    }
    sign = 0;
    if (a<0) {
      sign = 1;
      a = -a; // Also works for 0x80000000
    }
    // Normalize int
    int t=a; t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
    t = clz_tab[(t*clz_magic)>>>27]-1;

    exponent = 0x4000001E-t;
    mantissa = ((long)a)<<(32+t);
    return this;
  }

  /**
   * Assigns this <code>Real</code> the value of a signed long integer.
   * All long values can be represented without loss of accuracy.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = (double)a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>long</code> to assign.
   */
  public void assign(long a) {
    sign = 0;
    if (a<0) {
      sign = 1;
      a = -a; // Also works for 0x8000000000000000
    }
    exponent = 0x4000003E;
    mantissa = a;
    normalize();
  }

  /**
   * Assigns this <code>Real</code> a value encoded in a <code>String</code>
   * using base-10, as specified in {@link #assign(String,int)}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Double.{@link Double#valueOf(String) valueOf}(a);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     �-1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     100
   * </td></tr></table>
   *
   * @param a the <code>String</code> to assign.
   */
  public void assign(String a) {
    assign(a,10);
  }

  /**
   * Assigns this <code>Real</code> a value encoded in a <code>String</code>
   * using the specified number base. The string is parsed as follows:
   *
   * <ul>
   *   <li>If the string is <code>null</code> or an empty string, zero is
   *       assigned.
   *   <li>Leading spaces are ignored.
   *   <li>An optional sign, '+', '-' or '/', where '/' precedes a negative
   *       two's-complement number, reading: "an infinite number of 1-bits
   *       preceding the number".
   *   <li>Optional digits preceding the radix, in the specified base.
   *       <ul>
   *           <li>In base-2, allowed digits are '01'.
   *           <li>In base-8, allowed digits are '01234567'.
   *           <li>In base-10, allowed digits are '0123456789'.
   *           <li>In base-16, allowed digits are '0123456789ABCDEF'.
           </ul>
   *   <li>An optional radix character, '.' or ','.
   *   <li>Optional digits following the radix.
   *   <li>The following spaces are ignored.
   *   <li>An optional exponent indicator, 'e'. If not base-16, or after a
   *       space, 'E' is also accepted.
   *   <li>An optional sign, '+' or '-'.
   *   <li>Optional exponent digits <i><b>in base-10</b></i>.
   * </ul>
   *
   * <p><i>Valid examples:</i><br>
   * &nbsp;&nbsp;&nbsp;&nbsp;base-2:  <code>"-.110010101e+5"</code><br>
   * &nbsp;&nbsp;&nbsp;&nbsp;base-8:  <code>"+5462E-99"</code><br>
   * &nbsp;&nbsp;&nbsp;&nbsp;base-10: <code>"&nbsp;&nbsp;3,1415927"</code><br>
   * &nbsp;&nbsp;&nbsp;&nbsp;base-16: <code>"/FFF800C.CCCE e64"</code>
   *
   * <p>The number is parsed until the end of the string or an unknown
   * character is encountered, then silently returns even if the whole
   * string has not been parsed. Please note that specifying an
   * excessive number of digits in base-10 may in fact decrease the
   * accuracy of the result because of the extra multiplications performed.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td>
   * <td colspan="2">
   *     <code>this = Double.{@link Double#valueOf(String) valueOf}(a);
   *           // Works only for base-10</code>
   * </td></tr><tr><td valign="top" rowspan="2"><i>
   * Approximate&nbsp;error&nbsp;bound:</i>
   * </td><td width="1%">base-10</td><td>
   *     �-1 ULPs
   * </td></tr><tr><td>2/8/16</td><td>
   *     � ULPs
   * </td></tr><tr><td valign="top" rowspan="4"><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;</i>
   * </td><td width="1%">base-2</td><td>
   *     54
   * </td></tr><tr><td>base-8</td><td>
   *     60
   * </td></tr><tr><td>base-10</td><td>
   *     100
   * </td></tr><tr><td>base-16&nbsp;&nbsp;</td><td>
   *     60
   * </td></tr></table>
   *
   * @param a the <code>String</code> to assign.
   * @param base the number base of <code>a</code>. Valid base values are
   *     2, 8, 10 and 16.
   */
  public void assign(String a, int base) {
    if (a==null || a.length()==0) {
      assign(ZERO);
      return;
    }
    atof(a,base);
  }

  /**
   * Assigns this <code>Real</code> a value by directly setting the
   * fields of the internal representation. The arguments must
   * represent a valid, normalized <code>Real</code>. This is the
   * fastest way of assigning a constant value.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = (1-2*s) * m *
   *           Math.{@link Math#pow(double,double) pow}(2.0,e-62);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @param s {@link #sign} bit, 0 for positive sign, 1 for negative sign
   * @param e {@link #exponent}
   * @param m {@link #mantissa}
   */
  public void assign(int s, int e, long m) {
    sign = (byte)s;
    exponent = e;
    mantissa = m;
  }

  /**
   * Assigns this <code>Real</code> a value previously encoded into
   * into twelve consecutive bytes in a byte array using
   * {@link #toBytes(byte[],int) toBytes}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.2
   * </td></tr></table>
   *
   * @param data byte array to decode into this <code>Real</code>.
   * @param offset offset to start encoding from. The bytes
   *     <code>data[offset]...data[offset+11]</code> will be
   *     read.
   */
  public void assign(byte [] data, int offset) {
    sign = (byte)((data[offset+4]>>7)&1);
    exponent = (((data[offset ]&0xff)<<24)+
                ((data[offset +1]&0xff)<<16)+
                ((data[offset +2]&0xff)<<8)+
                ((data[offset +3]&0xff)));
    mantissa = (((long)(data[offset+ 4]&0x7f)<<56)+
                ((long)(data[offset+ 5]&0xff)<<48)+
                ((long)(data[offset+ 6]&0xff)<<40)+
                ((long)(data[offset+ 7]&0xff)<<32)+
                ((long)(data[offset+ 8]&0xff)<<24)+
                ((long)(data[offset+ 9]&0xff)<<16)+
                ((long)(data[offset+10]&0xff)<< 8)+
                ((long)(data[offset+11]&0xff)));
  }

  /**
   * Encodes an accurate representation of this <code>Real</code> value
   * into twelve consecutive bytes in a byte array. Can be decoded using
   * {@link #assign(byte[],int)}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.2
   * </td></tr></table>
   *
   * @param data byte array to save this <code>Real</code> in.
   * @param offset offset to start encoding to. The bytes
   *     <code>data[offset]...data[offset+11]</code> will be
   *     written.
   */
  public void toBytes(byte [] data, int offset) {
    data[offset ] = (byte)(exponent>>24);
    data[offset+ 1] = (byte)(exponent>>16);
    data[offset+ 2] = (byte)(exponent>>8);
    data[offset+ 3] = (byte)(exponent);
    data[offset+ 4] = (byte)((sign<<7)+(mantissa>>56));
    data[offset+ 5] = (byte)(mantissa>>48);
    data[offset+ 6] = (byte)(mantissa>>40);
    data[offset+ 7] = (byte)(mantissa>>32);
    data[offset+ 8] = (byte)(mantissa>>24);
    data[offset+ 9] = (byte)(mantissa>>16);
    data[offset+10] = (byte)(mantissa>>8);
    data[offset+11] = (byte)(mantissa);
  }

  /**
   * Assigns this <code>Real</code> the value corresponding to a given
   * bit representation. The argument is considered to be a
   * representation of a floating-point value according to the IEEE
   * 754 floating-point "single format" bit layout.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>float</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Float.{@link Float#intBitsToFloat(int)
   *           intBitsToFloat}(bits);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.6
   * </td></tr></table>
   *
   * @param bits a <code>float</code> value encoded in an <code>int</code>.
   */
  public void assignFloatBits(int bits) {
    sign = (byte)(bits>>>31);
    exponent = (bits>>23)&0xff;
    mantissa = (long)(bits&0x007fffff)<<39;
    if (exponent == 0 && mantissa == 0)
      return; // Valid zero
    if (exponent == 0 && mantissa != 0) {
      // Degenerate small float
      exponent = 0x40000000-126;
      normalize();
      return;
    }
    if (exponent <= 254) {
      // Normal IEEE 754 float
      exponent += 0x40000000-127;
      mantissa |= 1L<<62;
      return;
    }
    if (mantissa == 0)
      makeInfinity(sign);
    else
      makeNan();
  }

  /**
   * Assigns this <code>Real</code> the value corresponding to a given
   * bit representation. The argument is considered to be a
   * representation of a floating-point value according to the IEEE
   * 754 floating-point "double format" bit layout.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Double.{@link Double#longBitsToDouble(long)
   *           longBitsToDouble}(bits);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.6
   * </td></tr></table>
   *
   * @param bits a <code>double</code> value encoded in a <code>long</code>.
   */
  public void assignDoubleBits(long bits) {
    sign = (byte)((bits>>63)&1);
    exponent = (int)((bits>>52)&0x7ff);
    mantissa = (bits&0x000fffffffffffffL)<<10;
    if (exponent == 0 && mantissa == 0)
      return; // Valid zero
    if (exponent == 0 && mantissa != 0) {
      // Degenerate small float
      exponent = 0x40000000-1022;
      normalize();
      return;
    }
    if (exponent <= 2046) {
      // Normal IEEE 754 float
      exponent += 0x40000000-1023;
      mantissa |= 1L<<62;
      return;
    }
    if (mantissa == 0)
      makeInfinity(sign);
    else
      makeNan();
  }

  /**
   * Returns a representation of this <code>Real</code> according to
   * the IEEE 754 floating-point "single format" bit layout.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>float</code><i>&nbsp;code:</i></td><td>
   *     <code>Float.{@link Float#floatToIntBits(float)
   *           floatToIntBits}(this)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.7
   * </td></tr></table>
   *
   * @return the bits that represent the floating-point number.
   */
  public int toFloatBits() {
    if ((this.exponent < 0 && this.mantissa != 0))
      return 0x7fffffff; // nan
    int e = exponent-0x40000000+127;
    long m = mantissa;
    // Round properly!
    m += 1L<<38;
    if (m<0) {
      m >>>= 1;
      e++;
      if (exponent < 0) // Overflow
        return (sign<<31)|0x7f800000; // inf
    }
    if ((this.exponent < 0 && this.mantissa == 0) || e > 254)
      return (sign<<31)|0x7f800000; // inf
    if ((this.exponent == 0 && this.mantissa == 0) || e < -22)
      return (sign<<31); // zero
    if (e > 0)
      // Normal IEEE 754 float
      return (sign<<31)|(e<<23)|((int)(m>>>39)&0x007fffff);
    // Degenerate small float
    return (sign<<31)|((int)(m>>>(40-e))&0x007fffff);
  }

  /**
   * Returns a representation of this <code>Real</code> according to
   * the IEEE 754 floating-point "double format" bit layout.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>Double.{@link Double#doubleToLongBits(double)
   *           doubleToLongBits}(this)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.7
   * </td></tr></table>
   *
   * @return the bits that represent the floating-point number.
   */
  public long toDoubleBits() {
    if ((this.exponent < 0 && this.mantissa != 0))
      return 0x7fffffffffffffffL; // nan
    int e = exponent-0x40000000+1023;
    long m = mantissa;
    // Round properly!
    m += 1L<<9;
    if (m<0) {
      m >>>= 1;
      e++;
      if (exponent < 0)
        return ((long)sign<<63)|0x7ff0000000000000L; // inf
    }
    if ((this.exponent < 0 && this.mantissa == 0) || e > 2046)
      return ((long)sign<<63)|0x7ff0000000000000L; // inf
    if ((this.exponent == 0 && this.mantissa == 0) || e < -51)
      return ((long)sign<<63); // zero
    if (e > 0)
      // Normal IEEE 754 double
      return ((long)sign<<63)|((long)e<<52)|((m>>>10)&0x000fffffffffffffL);
    // Degenerate small double
    return ((long)sign<<63)|((m>>>(11-e))&0x000fffffffffffffL);
  }

  /**
   * Makes this <code>Real</code> the value of positive zero.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = 0;</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.2
   * </td></tr></table>
   */
  public void makeZero() {
    sign = 0;
    mantissa = 0;
    exponent = 0;
  }

  /**
   * Makes this <code>Real</code> the value of zero with the specified
   * sign.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = 0.0 * (1-2*s);</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.2
   * </td></tr></table>
   *
   * @param s sign bit, 0 to make a positive zero, 1 to make a
   * negative zero
   */
  public void makeZero(int s) {
    sign = (byte)s;
    mantissa = 0;
    exponent = 0;
  }

  /**
   * Makes this <code>Real</code> the value of infinity with the
   * specified sign.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Double.{@link Double#POSITIVE_INFINITY POSITIVE_INFINITY}
   *           * (1-2*s);</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @param s sign bit, 0 to make positive infinity, 1 to make
   * negative infinity
   */
  public void makeInfinity(int s) {
    sign = (byte)s;
    mantissa = 0;
    exponent = 0x80000000;
  }

  /**
   * Makes this <code>Real</code> the value of Not-a-Number (NaN).
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Double.{@link Double#NaN NaN};</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   */
  public void makeNan() {
    sign = 0;
    mantissa = 0x4000000000000000L;
    exponent = 0x80000000;
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * is zero, <code>false</code> otherwise.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this == 0)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object is
   *     zero, <code>false</code> otherwise.
   */
  public boolean isZero() {
    return (exponent == 0 && mantissa == 0);
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * is infinite, <code>false</code> otherwise.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>Double.{@link Double#isInfinite(double) isInfinite}(this)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object is
   *     infinite, <code>false</code> if it is finite or NaN.
   */
  public boolean isInfinity() {
    return (exponent < 0 && mantissa == 0);
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * is Not-a-Number (NaN), <code>false</code> otherwise.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>Double.{@link Double#isNaN(double) isNaN}(this)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object is
   *     NaN, <code>false</code> otherwise.
   */
  public boolean isNan() {
    return (exponent < 0 && mantissa != 0);
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * is finite, <code>false</code> otherwise.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(!Double.{@link Double#isNaN(double) isNaN}(this) &&
   *            !Double.{@link Double#isInfinite(double)
   *            isInfinite}(this))</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object is
   *     finite, <code>false</code> if it is infinite or NaN.
   */
  public boolean isFinite() {
    // That is, non-infinite and non-nan
    return (exponent >= 0);
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * is finite and nonzero, <code>false</code> otherwise.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *  <code>(!Double.{@link Double#isNaN(double) isNaN}(this) &&
   *         !Double.{@link Double#isInfinite(double) isInfinite}(this) &&
   *         (this!=0))</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object is
   *     finite and nonzero, <code>false</code> if it is infinite, NaN or zero.
   */
  public boolean isFiniteNonZero() {
    // That is, non-infinite and non-nan and non-zero
    return (exponent >= 0 && mantissa != 0);
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * is negative, <code>false</code> otherwise.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this < 0)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object
   *     is negative, <code>false</code> if it is positive or NaN.
   */
  public boolean isNegative() {
    return sign!=0;
  }

  /**
   * Calculates the absolute value.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#abs(double) abs}(this);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.2
   * </td></tr></table>
   */
  public void abs() {
    sign = 0;
  }

  /**
   * Negates this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = -this;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.2
   * </td></tr></table>
   */
  public void neg() {
    if (!(this.exponent < 0 && this.mantissa != 0))
      sign ^= 1;
  }

  /**
   * Copies the sign from <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#abs(double)
   *           abs}(this)*Math.{@link Math#signum(double) signum}(a);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.2
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to copy the sign from.
   */
  public void copysign(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    sign = a.sign;
  }

  /**
   * Readjusts the mantissa of this <code>Real</code>. The exponent is
   * adjusted accordingly. This is necessary when the mantissa has
   * been {@link #mantissa modified directly} for some purpose and may
   * be denormalized. The normalized mantissa of a finite
   * <code>Real</code> must have bit 63 cleared and bit 62 set. Using
   * a denormalized <code>Real</code> in <u>any</u> other operation
   * may produce undefined results.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.7
   * </td></tr></table>
   */
  public void normalize() {
    if ((this.exponent >= 0)) {
      if (mantissa > 0)
      {
        int clz = 0;
        int t = (int)(mantissa>>>32);
        if (t == 0) { clz = 32; t = (int)mantissa; }
        t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
        clz += clz_tab[(t*clz_magic)>>>27]-1;

        mantissa <<= clz;
        exponent -= clz;
        if (exponent < 0) // Underflow
          makeZero(sign);
      }
      else if (mantissa < 0)
      {
        mantissa = (mantissa+1)>>>1;
        exponent ++;
        if (mantissa == 0) { // Ooops, it was 0xffffffffffffffffL
          mantissa = 0x4000000000000000L;
          exponent ++;
        }
        if (exponent < 0) // Overflow
          makeInfinity(sign);
      }
      else // mantissa == 0
      {
        exponent = 0;
      }
    }
  }

  /**
   * Readjusts the mantissa of a <code>Real</code> with extended
   * precision. The exponent is adjusted accordingly. This is
   * necessary when the mantissa has been {@link #mantissa modified
   * directly} for some purpose and may be denormalized. The
   * normalized mantissa of a finite <code>Real</code> must have bit
   * 63 cleared and bit 62 set. Using a denormalized <code>Real</code>
   * in <u>any</u> other operation may produce undefined results.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2<sup>-64</sup> ULPs (i.e. of a normal precision <code>Real</code>)
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.7
   * </td></tr></table>
   *
   * @param extra the extra 64 bits of mantissa of this extended precision
   *     <code>Real</code>.
   * @return the extra 64 bits of mantissa of the resulting extended precision
   *     <code>Real</code>.
   */
  public long normalize128(long extra) {
    if (!(this.exponent >= 0))
      return 0;
    if (mantissa == 0) {
      if (extra == 0) {
        exponent = 0;
        return 0;
      }
      mantissa = extra;
      extra = 0;
      exponent -= 64;
      if (exponent < 0) { // Underflow
        makeZero(sign);
        return 0;
      }
    }
    if (mantissa < 0) {
      extra = (mantissa<<63)+(extra>>>1);
      mantissa >>>= 1;
      exponent ++;
      if (exponent < 0) { // Overflow
        makeInfinity(sign);
        return 0;
      }
      return extra;
    }

    int clz = 0;
    int t = (int)(mantissa>>>32);
    if (t == 0) { clz = 32; t = (int)mantissa; }
    t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
    clz += clz_tab[(t*clz_magic)>>>27]-1;
    if (clz == 0)
      return extra;

    mantissa = (mantissa<<clz)+(extra>>>(64-clz));
    extra <<= clz;
    exponent -= clz;
    if (exponent < 0) { // Underflow
      makeZero(sign);
      return 0;
    }
    return extra;
  }

  /**
   * Rounds an extended precision <code>Real</code> to the nearest
   * <code>Real</code> of normal precision.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param extra the extra 64 bits of mantissa of this extended precision
   *     <code>Real</code>.
   */
  public void roundFrom128(long extra) {
    mantissa += (extra>>63)&1;
    normalize();
  }

  private int compare(Real a) {
    // Compare of normal floats, zeros, but not nan or equal-signed inf
    if ((this.exponent == 0 && this.mantissa == 0) && (a.exponent == 0 && a.mantissa == 0))
      return 0;
    if (sign != a.sign)
      return a.sign-sign;
    int s = (this.sign==0) ? 1 : -1;
    if ((this.exponent < 0 && this.mantissa == 0))
      return s;
    if ((a.exponent < 0 && a.mantissa == 0))
      return -s;
    if (exponent != a.exponent)
      return exponent<a.exponent ? -s : s;
    if (mantissa != a.mantissa)
      return mantissa<a.mantissa ? -s : s;
    return 0;
  }

  private boolean invalidCompare(Real a) {
    return ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0) ||
            ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign == a.sign));
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is equal to
   * <code>a</code>.
   * If the numbers are incomparable, i.e. the values are infinities of
   * the same sign or any of them is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this == a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is
   *     equal to the value represented by <code>a</code>. <code>false</code>
   *     otherwise, or if the numbers are incomparable.
   */
  public boolean equalTo(Real a) {
    if (invalidCompare(a))
      return false;
    return compare(a) == 0;
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is not equal to
   * <code>a</code>.
   * If the numbers are incomparable, i.e. the values are infinities of
   * the same sign or any of them is NaN, <code>false</code> is always
   * returned.
   * This distinguishes <code>notEqualTo(a)</code> from the expression
   * <code>!equalTo(a)</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this != a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is not
   *     equal to the value represented by <code>a</code>. <code>false</code>
   *     otherwise, or if the numbers are incomparable.
   */
  public boolean notEqualTo(Real a) {
    if (invalidCompare(a))
      return false;
    return compare(a) != 0;
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is less than
   * <code>a</code>.
   * If the numbers are incomparable, i.e. the values are infinities of
   * the same sign or any of them is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &lt; a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is less
   *     than the value represented by <code>a</code>. <code>false</code>
   *     otherwise, or if the numbers are incomparable.
   */
  public boolean lessThan(Real a) {
    if (invalidCompare(a))
      return false;
    return compare(a) < 0;
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is less than or
   * equal to <code>a</code>.
   * If the numbers are incomparable, i.e. the values are infinities of
   * the same sign or any of them is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &lt;= a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is less
   *     than or equal to the value represented by <code>a</code>.
   *     <code>false</code> otherwise, or if the numbers are incomparable.
   */
  public boolean lessEqual(Real a) {
    if (invalidCompare(a))
      return false;
    return compare(a) <= 0;
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is greater than
   * <code>a</code>.
   * If the numbers are incomparable, i.e. the values are infinities of
   * the same sign or any of them is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &gt; a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is
   *     greater than the value represented by <code>a</code>.
   *     <code>false</code> otherwise, or if the numbers are incomparable.
   */
  public boolean greaterThan(Real a) {
    if (invalidCompare(a))
      return false;
    return compare(a) > 0;
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is greater than
   * or equal to <code>a</code>.
   * If the numbers are incomparable, i.e. the values are infinities of
   * the same sign or any of them is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &gt;= a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is
   *     greater than or equal to the value represented by <code>a</code>.
   *     <code>false</code> otherwise, or if the numbers are incomparable.
   */
  public boolean greaterEqual(Real a) {
    if (invalidCompare(a))
      return false;
    return compare(a) >= 0;
  }

  /**
   * Returns <code>true</code> if the absolute value of this
   * <code>Real</code> is less than the absolute value of
   * <code>a</code>.
   * If the numbers are incomparable, i.e. the values are both infinite
   * or any of them is NaN, <code>false</code> is always returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(Math.{@link Math#abs(double) abs}(this) &lt;
   *           Math.{@link Math#abs(double) abs}(a))</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.5
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to compare to this.
   * @return <code>true</code> if the absolute of the value represented by
   *     this object is less  than the absolute of the value represented by
   *     <code>a</code>.
   *     <code>false</code> otherwise, or if the numbers are incomparable.
   */
  public boolean absLessThan(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0) || (this.exponent < 0 && this.mantissa == 0))
      return false;
    if ((a.exponent < 0 && a.mantissa == 0))
      return true;
    if (exponent != a.exponent)
      return exponent<a.exponent;
    return mantissa<a.mantissa;
  }

  /**
   * Multiplies this <code>Real</code> by 2 to the power of <code>n</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * This operation is faster than normal multiplication since it only
   * involves adding to the exponent.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this *= Math.{@link Math#pow(double,double) pow}(2.0,n);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.3
   * </td></tr></table>
   */
  public void scalbn(int n) {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    exponent += n;
    if (exponent < 0) {
      if (n<0)
        makeZero(sign); // Underflow
      else
        makeInfinity(sign); // Overflow
    }
  }

  /**
   * Calculates the next representable neighbour of this <code>Real</code>
   * in the direction towards <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * If the two values are equal, nothing happens.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this += Math.{@link Math#ulp(double) ulp}(this)*Math.{@link
   *           Math#signum(double) signum}(a-this);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.8
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument.
   */
  public void nextafter(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign == a.sign)
      return;
    int dir = -compare(a);
    if (dir == 0)
      return;
    if ((this.exponent == 0 && this.mantissa == 0)) {
      { this.mantissa = MIN.mantissa; this.exponent = MIN.exponent; this.sign = MIN.sign; };
      sign = (byte)(dir<0 ? 1 : 0);
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0)) {
      { this.mantissa = MAX.mantissa; this.exponent = MAX.exponent; this.sign = MAX.sign; };
      sign = (byte)(dir<0 ? 0 : 1);
      return;
    }
    if ((this.sign==0) ^ dir<0) {
      mantissa ++;
    } else {
      if (mantissa == 0x4000000000000000L) {
        mantissa <<= 1;
        exponent--;
      }
      mantissa --;
    }
    normalize();
  }

  /**
   * Calculates the largest (closest to positive infinity)
   * <code>Real</code> value that is less than or equal to this
   * <code>Real</code> and is equal to a mathematical integer.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#floor(double) floor}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.5
   * </td></tr></table>
   */
  public void floor() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    if (exponent < 0x40000000) {
      if ((this.sign==0))
        makeZero(sign);
      else {
        exponent = ONE.exponent;
        mantissa = ONE.mantissa;
        // sign unchanged!
      }
      return;
    }
    int shift = 0x4000003e-exponent;
    if (shift <= 0)
      return;

    if ((this.sign!=0))
      mantissa += ((1L<<shift)-1);

    mantissa &= ~((1L<<shift)-1);

    if ((this.sign!=0))
      normalize();
  }

  /**
   * Calculates the smallest (closest to negative infinity)
   * <code>Real</code> value that is greater than or equal to this
   * <code>Real</code> and is equal to a mathematical integer.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#ceil(double) ceil}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.8
   * </td></tr></table>
   */
  public void ceil() {
    neg();
    floor();
    neg();
  }

  /**
   * Rounds this <code>Real</code> value to the closest value that
   * is equal to a mathematical integer. If two
   * <code>Real</code> values that are mathematical integers are equally close,
   * the result is the integer value with the largest magnitude (positive or
   * negative).
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#rint(double) rint}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.3
   * </td></tr></table>
   */
  public void round() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    if (exponent < 0x3fffffff) {
      makeZero(sign);
      return;
    }
    int shift = 0x4000003e-exponent;
    if (shift <= 0)
      return;

    mantissa += 1L<<(shift-1); // Bla-bla, this works almost
    mantissa &= ~((1L<<shift)-1);
    normalize();
  }

  /**
   * Rounds this <code>Real</code> value to the closest value towards
   * zero that is equal to a mathematical integer.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = (double)((long)this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.2
   * </td></tr></table>
   */
  public void trunc() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    if (exponent < 0x40000000) {
      makeZero(sign);
      return;
    }
    int shift = 0x4000003e-exponent;
    if (shift <= 0)
      return;

    mantissa &= ~((1L<<shift)-1);
    normalize();
  }

  /**
   * Calculates the fractional part of this <code>Real</code> by subtracting
   * the closest value towards zero that is equal to a mathematical integer.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this -= (double)((long)this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.2
   * </td></tr></table>
   */
  public void frac() {
    if (!(this.exponent >= 0 && this.mantissa != 0) || exponent < 0x40000000)
      return;

    int shift = 0x4000003e-exponent;
    if (shift <= 0) {
      makeZero(sign);
      return;
    }

    mantissa &= ((1L<<shift)-1);
    normalize();
  }

  /**
   * Converts this <code>Real</code> value to the closest
   * <code>int</code> value towards zero.
   *
   * <p>If the value of this <code>Real</code> is too large, {@link
   * Integer#MAX_VALUE} is returned. However, if the value of this
   * <code>Real</code> is too small, <code>-Integer.MAX_VALUE</code>
   * is returned, not {@link Integer#MIN_VALUE}. This is done to
   * ensure that the sign will be correct if you calculate
   * <code>-this.toInteger()</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(int)this</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.6
   * </td></tr></table>
   *
   * @return an <code>int</code> representation of this <code>Real</code>.
   */
  public int toInteger() {
    if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa != 0))
      return 0;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      if ((this.sign==0))
        return 0x7fffffff;
      else
        return 0x80000001; // So that you can take -x.toInteger()
    }
    if (exponent < 0x40000000)
      return 0;
    int shift = 0x4000003e-exponent;
    if (shift < 32) {
      if ((this.sign==0))
        return 0x7fffffff;
      else
        return 0x80000001; // So that you can take -x.toInteger()
    }
    return (this.sign==0) ?
      (int)(mantissa>>>shift) : -(int)(mantissa>>>shift);
  }

  /**
   * Converts this <code>Real</code> value to the closest
   * <code>long</code> value towards zero.
   *
   * <p>If the value of this <code>Real</code> is too large, {@link
   * Long#MAX_VALUE} is returned. However, if the value of this
   * <code>Real</code> is too small, <code>-Long.MAX_VALUE</code> is
   * returned, not {@link Long#MIN_VALUE}. This is done to ensure that
   * the sign will be correct if you calculate
   * <code>-this.toLong()</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(long)this</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.5
   * </td></tr></table>
   *
   * @return a <code>long</code> representation of this <code>Real</code>.
   */
  public long toLong() {
    if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa != 0))
      return 0;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      if ((this.sign==0))
        return 0x7fffffffffffffffL;
      else
        return 0x8000000000000001L; // So that you can take -x.toLong()
    }
    if (exponent < 0x40000000)
      return 0;
    int shift = 0x4000003e-exponent;
    if (shift < 0) {
      if ((this.sign==0))
        return 0x7fffffffffffffffL;
      else
        return 0x8000000000000001L; // So that you can take -x.toLong()
    }
    return (this.sign==0) ? (mantissa>>>shift) : -(mantissa>>>shift);
  }

  /**
   * Returns <code>true</code> if the value of this <code>Real</code>
   * represents a mathematical integer. If the value is too large to
   * determine if it is an integer, <code>true</code> is returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this == (long)this)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.6
   * </td></tr></table>
   *
   * @return <code>true</code> if the value represented by this object
   *     represents a mathematical integer, <code>false</code> otherwise.
   */
  public boolean isIntegral() {
    if ((this.exponent < 0 && this.mantissa != 0))
      return false;
    if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa == 0))
      return true;
    if (exponent < 0x40000000)
      return false;
    int shift = 0x4000003e-exponent;
    if (shift <= 0)
      return true;
    return (mantissa&((1L<<shift)-1)) == 0;
  }

  /**
   * Returns <code>true</code> if the mathematical integer represented
   * by this <code>Real</code> is odd. You <u>must</u> first determine
   * that the value is actually an integer using {@link
   * #isIntegral()}. If the value is too large to determine if the
   * integer is odd, <code>false</code> is returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>((((long)this)&1) == 1)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.6
   * </td></tr></table>
   *
   * @return <code>true</code> if the mathematical integer represented by
   *     this <code>Real</code> is odd, <code>false</code> otherwise.
   */
  public boolean isOdd() {
    if (!(this.exponent >= 0 && this.mantissa != 0) ||
        exponent < 0x40000000 || exponent > 0x4000003e)
      return false;
    int shift = 0x4000003e-exponent;
    return ((mantissa>>>shift)&1) != 0;
  }

  /**
   * Exchanges the contents of this <code>Real</code> and <code>a</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>tmp=this; this=a; a=tmp;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     0.5
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to exchange with this.
   */
  public void swap(Real a) {
    long tmpMantissa=mantissa; mantissa=a.mantissa; a.mantissa=tmpMantissa;
    int tmpExponent=exponent; exponent=a.exponent; a.exponent=tmpExponent;
    byte tmpSign =sign; sign =a.sign; a.sign =tmpSign;
  }

  // Temporary values used by functions (to avoid "new" inside functions)
  private static Real tmp0 = new Real(); // tmp for basic functions
  private static Real recipTmp = new Real();
  private static Real recipTmp2 = new Real();
  private static Real sqrtTmp = new Real();
  private static Real expTmp = new Real();
  private static Real expTmp2 = new Real();
  private static Real expTmp3 = new Real();
  private static Real tmp1 = new Real();
  private static Real tmp2 = new Real();
  private static Real tmp3 = new Real();
  private static Real tmp4 = new Real();
  private static Real tmp5 = new Real();

  /**
   * Calculates the sum of this <code>Real</code> and <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this += a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     �� 1.0 ��
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to add to this.
   */
  public Real add(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return this;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign != a.sign)
        makeNan();
      else
        makeInfinity((this.exponent < 0 && this.mantissa == 0) ? sign : a.sign);
      return this;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      if ((this.exponent == 0 && this.mantissa == 0))
        { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
      if ((this.exponent == 0 && this.mantissa == 0))
        sign=0;
      return this;
    }

    byte s;
    int e;
    long m;
    if (exponent > a.exponent ||
        (exponent == a.exponent && mantissa>=a.mantissa))
    {
      s = a.sign;
      e = a.exponent;
      m = a.mantissa;
    } else {
      s = sign;
      e = exponent;
      m = mantissa;
      sign = a.sign;
      exponent = a.exponent;
      mantissa = a.mantissa;
    }
    int shift = exponent-e;
    if (shift>=64)
      return this;

    if (sign == s) {
      mantissa += m>>>shift;
      if (mantissa >= 0 && shift>0 && ((m>>>(shift-1))&1) != 0)
        mantissa ++; // We don't need normalization, so round now
      if (mantissa < 0) {
        // Simplified normalize()
        mantissa = (mantissa+1)>>>1;
        exponent ++;
        if (exponent < 0) { // Overflow
          makeInfinity(sign);
          return this;
        }
      }
    } else {
      if (shift>0) {
        // Shift mantissa up to increase accuracy
        mantissa <<= 1;
        exponent --;
        shift --;
      }
      m = -m;
      mantissa += m>>shift;
      if (mantissa >= 0 && shift>0 && ((m>>>(shift-1))&1) != 0)
        mantissa ++; // We don't need to shift down, so round now
      if (mantissa < 0) {
        // Simplified normalize()
        mantissa = (mantissa+1)>>>1;
        exponent ++; // Can't overflow
      } else if (shift==0) {
        // Operands have similar magnitude => many bits may be cancelled

        // Magic rounding: if result of subtract leaves only a few bits
        // standing, the result should most likely be 0...
        if (magicRounding && mantissa > 0 && mantissa <= 7)
          mantissa = 0;
        normalize();
      } // else... if (shift>=1 && mantissa>=0) it should be a-ok
    }

    if ((this.exponent == 0 && this.mantissa == 0))
      sign=0;
    
    return this;
  }

  /**
   * Calculates the sum of this <code>Real</code> and <code>a</code> with
   * extended precision.
   * Replaces the contents of this <code>Real</code> with the result.
   * Returns the extra mantissa of the extended precision result.
   *
   * <p>An extra 64 bits of mantissa is added to both arguments for
   * extended precision. If any of the arguments are not of extended
   * precision, use <code>0</code> for the extra mantissa.
   *
   * <p>Extended prevision can be useful in many situations. For
   * instance, when accumulating a lot of very small values it is
   * advantageous for the accumulator to have extended precision. To
   * convert the extended precision value back to a normal
   * <code>Real</code> for further processing, use {@link
   * #roundFrom128(long)}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this += a;</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2<sup>-62</sup> ULPs (i.e. of a normal precision <code>Real</code>)
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     2.0
   * </td></tr></table>
   *
   * @param extra the extra 64 bits of mantissa of this extended precision
   *     <code>Real</code>.
   * @param a the <code>Real</code> to add to this.
   * @param aExtra the extra 64 bits of mantissa of the extended precision
   *     value <code>a</code>.
   * @return the extra 64 bits of mantissa of the resulting extended precision
   *     <code>Real</code>.
   */
  public long add128(long extra, Real a, long aExtra) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return 0;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign != a.sign)
        makeNan();
      else
        makeInfinity((this.exponent < 0 && this.mantissa == 0) ? sign : a.sign);
      return 0;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      if ((this.exponent == 0 && this.mantissa == 0)) {
        { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
        extra = aExtra;
      }
      if ((this.exponent == 0 && this.mantissa == 0))
        sign=0;
      return extra;
    }

    byte s;
    int e;
    long m;
    long x;
    if (exponent > a.exponent ||
        (exponent == a.exponent && mantissa>a.mantissa) ||
        (exponent == a.exponent && mantissa==a.mantissa &&
         (extra>>>1)>=(aExtra>>>1)))
    {
      s = a.sign;
      e = a.exponent;
      m = a.mantissa;
      x = aExtra;
    } else {
      s = sign;
      e = exponent;
      m = mantissa;
      x = extra;
      sign = a.sign;
      exponent = a.exponent;
      mantissa = a.mantissa;
      extra = aExtra;
    }

    int shift = exponent-e;
    if (shift>=127)
      return extra;
    if (shift>=64) {
      x = m>>>(shift-64);
      m = 0;
    } else if (shift>0) {
      x = (x>>>shift)+(m<<(64-shift));
      m >>>= shift;
    }
    extra >>>= 1;
    x >>>= 1;

    if (sign == s) {
      extra += x;
      mantissa += (extra>>63)&1;
      mantissa += m;
    } else {
      extra -= x;
      mantissa -= (extra>>63)&1;
      mantissa -= m;
      // Magic rounding: if result of subtract leaves only a few bits
      // standing, the result should most likely be 0...
      if (mantissa == 0 && extra > 0 && extra <= 0x1f)
        extra = 0;
    }
    extra <<= 1;

    extra = normalize128(extra);

    if ((this.exponent == 0 && this.mantissa == 0))
      sign=0;

    return extra;
  }

  /**
   * Calculates the difference between this <code>Real</code> and
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this -= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     2.0
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to subtract from this.
   */
  public void sub(Real a) {
    tmp0.mantissa = a.mantissa;
    tmp0.exponent = a.exponent;
    tmp0.sign = (byte)(a.sign^1);
    add(tmp0);
  }

  /**
   * Calculates the product of this <code>Real</code> and the integer
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this *= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.3
   * </td></tr></table>
   *
   * @param a the <code>int</code> to multiply to this.
   */
  public Real mul(int a) {
    if ((this.exponent < 0 && this.mantissa != 0))
      return this;
    if (a<0) {
      sign ^= 1;
      a = -a;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || a==0) {
      if ((this.exponent < 0 && this.mantissa == 0))
        makeNan();
      else
        makeZero(sign);
      return this;
    }
    if ((this.exponent < 0 && this.mantissa == 0))
      return this;

    // Normalize int
    int t=a; t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
    t = clz_tab[(t*clz_magic)>>>27];
    exponent += 0x1F-t;
    a <<= t;

    if (exponent < 0) {
      makeInfinity(sign); // Overflow
      return this;
    }

    long a0 = mantissa & 0x7fffffff;
    long a1 = mantissa >>> 31;
    long b0 = a & 0xffffffffL;

    mantissa = a1*b0;
    // If we're going to need normalization, we don't want to round twice
    int round = (mantissa<0) ? 0 : 0x40000000;
    mantissa += ((a0*b0 + round)>>>31);

    // Simplified normalize()
    if (mantissa < 0) {
      mantissa = (mantissa+1)>>>1;
      exponent ++;
      if (exponent < 0) // Overflow
        makeInfinity(sign);
    }
    return this;
  }

  /**
   * Calculates the product of this <code>Real</code> and <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this *= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.3
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to multiply to this.
   */
  public void mul(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    sign ^= a.sign;
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0))
        makeNan();
      else
        makeZero(sign);
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      makeInfinity(sign);
      return;
    }

    long a0 = mantissa & 0x7fffffff;
    long a1 = mantissa >>> 31;
    long b0 = a.mantissa & 0x7fffffff;
    long b1 = a.mantissa >>> 31;

    mantissa = a1*b1;
    // If we're going to need normalization, we don't want to round twice
    int round = (mantissa<0) ? 0 : 0x40000000;
    mantissa += ((a0*b1 + a1*b0 + ((a0*b0)>>>31) + round)>>>31);

    int aExp = a.exponent;
    exponent += aExp-0x40000000;
    if (exponent < 0) {
      if (exponent == -1 && aExp < 0x40000000 && mantissa < 0) {
        // Not underflow after all, it will be corrected in the
        // normalization below
      } else {
        if (aExp < 0x40000000)
          makeZero(sign); // Underflow
        else
          makeInfinity(sign); // Overflow
        return;
      }
    }

    // Simplified normalize()
    if (mantissa < 0) {
      mantissa = (mantissa+1)>>>1;
      exponent ++;
      if (exponent < 0) // Overflow
        makeInfinity(sign);
    }
  }

  /**
   * Calculates the product of this <code>Real</code> and <code>a</code> with
   * extended precision.
   * Replaces the contents of this <code>Real</code> with the result.
   * Returns the extra mantissa of the extended precision result.
   *
   * <p>An extra 64 bits of mantissa is added to both arguments for
   * extended precision. If any of the arguments are not of extended
   * precision, use <code>0</code> for the extra mantissa. See also {@link
   * #add128(long,Real,long)}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this *= a;</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2<sup>-60</sup> ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     3.1
   * </td></tr></table>
   *
   * @param extra the extra 64 bits of mantissa of this extended precision
   *     <code>Real</code>.
   * @param a the <code>Real</code> to multiply to this.
   * @param aExtra the extra 64 bits of mantissa of the extended precision
   *     value <code>a</code>.
   * @return the extra 64 bits of mantissa of the resulting extended precision
   *     <code>Real</code>.
   */
  public long mul128(long extra, Real a, long aExtra) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return 0;
    }
    sign ^= a.sign;
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0))
        makeNan();
      else
        makeZero(sign);
      return 0;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      makeInfinity(sign);
      return 0;
    }
    int aExp = a.exponent;
    exponent += aExp-0x40000000;
    if (exponent < 0) {
      if (aExp < 0x40000000)
        makeZero(sign); // Underflow
      else
        makeInfinity(sign); // Overflow
      return 0;
    }
    long ffffffffL = 0xffffffffL;
    long a0 = extra & ffffffffL;
    long a1 = extra >>> 32;
    long a2 = mantissa & ffffffffL;
    long a3 = mantissa >>> 32;
    long b0 = aExtra & ffffffffL;
    long b1 = aExtra >>> 32;
    long b2 = a.mantissa & ffffffffL;
    long b3 = a.mantissa >>> 32;

    a0 = ((a3*b0>>>2)+
          (a2*b1>>>2)+
          (a1*b2>>>2)+
          (a0*b3>>>2)+
          0x60000000)>>>28;
          //(a2*b0>>>34)+(a1*b1>>>34)+(a0*b2>>>34)+0x08000000)>>>28;
    a1 *= b3;
    b0 = a2*b2;
    b1 *= a3;

    a0 += ((a1<<2)&ffffffffL) + ((b0<<2)&ffffffffL) + ((b1<<2)&ffffffffL);
    a1 = (a0>>>32) + (a1>>>30) + (b0>>>30) + (b1>>>30);
    a0 &= ffffffffL;

    a2 *= b3;
    b2 *= a3;
    a1 += ((a2<<2)&ffffffffL) + ((b2<<2)&ffffffffL);

    extra = (a1<<32) + a0;
    mantissa = ((a3*b3)<<2) + (a1>>>32) + (a2>>>30) + (b2>>>30);

    extra = normalize128(extra);
    return extra;
  }

  private void mul10() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    mantissa += (mantissa+2)>>>2;
    exponent += 3;
    if (mantissa < 0) {
      mantissa = (mantissa+1)>>>1;
      exponent++;
    }
    if (exponent < 0)
      makeInfinity(sign); // Overflow
  }

  /**
   * Calculates the square of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = this*this;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.1
   * </td></tr></table>
   */
  public void sqr() {
    sign = 0;
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    int e = exponent;
    exponent += exponent-0x40000000;
    if (exponent < 0) {
      if (e < 0x40000000)
        makeZero(sign); // Underflow
      else
        makeInfinity(sign); // Overflow
      return;
    }
    long a0 = mantissa&0x7fffffff;
    long a1 = mantissa>>>31;

    mantissa = a1*a1;
    // If we're going to need normalization, we don't want to round twice
    int round = (mantissa<0) ? 0 : 0x40000000;
    mantissa += ((((a0*a1)<<1) + ((a0*a0)>>>31) + round)>>>31);

    // Simplified normalize()
    if (mantissa < 0) {
      mantissa = (mantissa+1)>>>1;
      exponent ++;
      if (exponent < 0) // Overflow
        makeInfinity(sign);
    }
  }

  /**
   * Calculates the quotient of this <code>Real</code> and the integer
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this /= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     2.6
   * </td></tr></table>
   *
   * @param a the <code>int</code> to divide this with.
   */
  public void div(int a) {
    if ((this.exponent < 0 && this.mantissa != 0))
      return;
    if (a<0) {
      sign ^= 1;
      a = -a;
    }
    if ((this.exponent < 0 && this.mantissa == 0))
      return;
    if ((this.exponent == 0 && this.mantissa == 0)) {
      if (a==0)
        makeNan();
      return;
    }
    if (a==0) {
      makeInfinity(sign);
      return;
    }
    long denom = a & 0xffffffffL;
    long remainder = mantissa%denom;
    mantissa /= denom;

    // Normalizing mantissa and scaling remainder accordingly
    int clz = 0;
    int t = (int)(mantissa>>>32);
    if (t == 0) { clz = 32; t = (int)mantissa; }
    t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
    clz += clz_tab[(t*clz_magic)>>>27]-1;
    mantissa <<= clz;
    remainder <<= clz;
    exponent -= clz;

    // Final division, correctly rounded
    remainder = (remainder+denom/2)/denom;
    mantissa += remainder;

    if (exponent < 0) // Underflow
      makeZero(sign);
  }

  private static long ldiv(long a, long b) {
    // Calculate (a<<63)/b, where a<2**64, b<2**63, b<=a and a<2*b
    // The result will always be 63 bits, leading to a 3-stage radix-2**21
    // (very high radix) algorithm, as described here:
    // S.F. Oberman and M.J. Flynn, "Division Algorithms and Implementations,"
    // IEEE Trans. Computers, vol. 46, no. 8, pp. 833-854, Aug 1997
    // Section 4: "Very High Radix Algorithms"

    int bInv24; // Approximate 1/b, never more than 24 bits
    int aHi24; // High 24 bits of a (sometimes 25 bits)
    int next21; // The next 21 bits of result, possibly 1 less
    long q; // Resulting quotient: round((a<<63)/b)
    // Preparations
    bInv24 = (int)(0x400000000000L/((b>>>40)+1));
    aHi24 = (int)(a>>32)>>>8;
    a <<= 20; // aHi24 and a overlap by 4 bits
    // Now perform the division
    next21 = (int)(((long)aHi24*(long)bInv24)>>>26);
    a -= next21*b; // Bits above 2**64 will always be cancelled
    // No need to remove remainder since this will be cared for in next block
    q = next21;
    aHi24 = (int)(a>>32)>>>7;
    a <<= 21;
    // Two more almost identical blocks...
    next21 = (int)(((long)aHi24*(long)bInv24)>>>26);
    a -= next21*b;
    q = (q<<21)+next21;
    aHi24 = (int)(a>>32)>>>7;
    a <<= 21;
    next21 = (int)(((long)aHi24*(long)bInv24)>>>26);
    a -= next21*b;
    q = (q<<21)+next21;
    // Remove final remainder
    if (a<0 || a>=b) { q++; a -= b; }
    a <<= 1;
    // Round correctly
    if (a<0 || a>=b) q++;
    return q;
  }

  /**
   * Calculates the quotient of this <code>Real</code> and <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this /= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     2.6
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to divide this with.
   */
  public void div(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    sign ^= a.sign;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      if ((a.exponent < 0 && a.mantissa == 0))
        makeNan();
      return;
    }
    if ((a.exponent < 0 && a.mantissa == 0)) {
      makeZero(sign);
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      if ((a.exponent == 0 && a.mantissa == 0))
        makeNan();
      return;
    }
    if ((a.exponent == 0 && a.mantissa == 0)) {
      makeInfinity(sign);
      return;
    }
    exponent += 0x40000000-a.exponent;
    if (mantissa < a.mantissa) {
      mantissa <<= 1;
      exponent--;
    }
    if (exponent < 0) {
      if (a.exponent >= 0x40000000)
        makeZero(sign); // Underflow
      else
        makeInfinity(sign); // Overflow
      return;
    }

    if (a.mantissa == 0x4000000000000000L)
      return;
    mantissa = ldiv(mantissa,a.mantissa);
  }

  /**
   * Calculates the quotient of <code>a</code> and this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = a/this;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     3.1
   * </td></tr></table>
   *
   * @param a the <code>Real</code> to be divided by this.
   */
  public void rdiv(Real a) {
    { recipTmp.mantissa = a.mantissa; recipTmp.exponent = a.exponent; recipTmp.sign = a.sign; };
    recipTmp.div(this);
    { this.mantissa = recipTmp.mantissa; this.exponent = recipTmp.exponent; this.sign = recipTmp.sign; };
  }

  /**
   * Calculates the reciprocal of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = 1/this;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     2.3
   * </td></tr></table>
   */
  public void recip() {
    if ((this.exponent < 0 && this.mantissa != 0))
      return;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      makeZero(sign);
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      makeInfinity(sign);
      return;
    }
    exponent = 0x80000000-exponent;
    if (mantissa == 0x4000000000000000L) {
      if (exponent < 0)
        makeInfinity(sign); // Overflow
      return;
    }
    exponent--;
    mantissa = ldiv(0x8000000000000000L,mantissa);
  }

  /**
   * Calculates the reciprocal of this <code>Real</code> with
   * extended precision.
   * Replaces the contents of this <code>Real</code> with the result.
   * Returns the extra mantissa of the extended precision result.
   *
   * <p>An extra 64 bits of mantissa is added for extended precision.
   * If the argument is not of extended precision, use <code>0</code>
   * for the extra mantissa. See also {@link #add128(long,Real,long)}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = 1/this;</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2<sup>-60</sup> ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     17
   * </td></tr></table>
   *
   * @param extra the extra 64 bits of mantissa of this extended precision
   *     <code>Real</code>.
   * @return the extra 64 bits of mantissa of the resulting extended precision
   *     <code>Real</code>.
   */
  public long recip128(long extra) {
    if ((this.exponent < 0 && this.mantissa != 0))
      return 0;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      makeZero(sign);
      return 0;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      makeInfinity(sign);
      return 0;
    }

    byte s = sign;
    sign = 0;

    // Special case, simple power of 2
    if (mantissa == 0x4000000000000000L && extra == 0) {
      exponent = 0x80000000-exponent;
      if (exponent<0) // Overflow
        makeInfinity(s);
      return 0;
    }

    // Normalize exponent
    int exp = 0x40000000-exponent;
    exponent = 0x40000000;

    // Save -A
    { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; };
    long recipTmpExtra = extra;
    recipTmp.neg();

    // First establish approximate result (actually 63 bit accurate)
    recip();

    // Perform one Newton-Raphson iteration
    // Xn+1 = Xn + Xn*(1-A*Xn)
    { recipTmp2.mantissa = this.mantissa; recipTmp2.exponent = this.exponent; recipTmp2.sign = this.sign; };
    extra = mul128(0,recipTmp,recipTmpExtra);
    extra = add128(extra,ONE,0);
    extra = mul128(extra,recipTmp2,0);
    extra = add128(extra,recipTmp2,0);

    // Fix exponent
    scalbn(exp);
    // Fix sign
    if (!isNan())
      sign = s;

    return extra;
  }

  /**
   * Calculates the mathematical integer that is less than or equal to
   * this <code>Real</code> divided by <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#floor(double) floor}(this/a);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     22
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument.
   */
  public void divf(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0)) {
      if ((a.exponent < 0 && a.mantissa == 0))
        makeNan();
      return;
    }
    if ((a.exponent < 0 && a.mantissa == 0)) {
      makeZero(sign);
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      if ((a.exponent == 0 && a.mantissa == 0))
        makeNan();
      return;
    }
    if ((a.exponent == 0 && a.mantissa == 0)) {
      makeInfinity(sign);
      return;
    }
    { tmp0.mantissa = a.mantissa; tmp0.exponent = a.exponent; tmp0.sign = a.sign; }; // tmp0 should be free
    // Perform same division as with mod, and don't round up
    long extra = tmp0.recip128(0);
    extra = mul128(0,tmp0,extra);
    if (((tmp0.sign!=0) && (extra < 0 || extra > 0x1f)) ||
        (!(tmp0.sign!=0) && extra < 0 && extra > 0xffffffe0))
    {
      // For accurate floor()
      mantissa++;
      normalize();
    }
    floor();
  }

  private void modInternal(/*long thisExtra,*/ Real a, long aExtra) {
    { tmp0.mantissa = a.mantissa; tmp0.exponent = a.exponent; tmp0.sign = a.sign; }; // tmp0 should be free
    long extra = tmp0.recip128(aExtra);
    extra = tmp0.mul128(extra,this,0/*thisExtra*/); // tmp0 == this/a
    if (tmp0.exponent > 0x4000003e) {
      // floor() will be inaccurate
      makeZero(a.sign); // What else can be done? makeNan?
      return;
    }
    if (((tmp0.sign!=0) && (extra < 0 || extra > 0x1f)) ||
        (!(tmp0.sign!=0) && extra < 0 && extra > 0xffffffe0))
    {
      // For accurate floor() with a bit of "magical rounding"
      tmp0.mantissa++;
      tmp0.normalize();
    }
    tmp0.floor();
    tmp0.neg(); // tmp0 == -floor(this/a)
    extra = tmp0.mul128(0,a,aExtra);
    extra = add128(0/*thisExtra*/,tmp0,extra);
    roundFrom128(extra);
  }

  /**
   * Calculates the value of this <code>Real</code> modulo <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The modulo in this case is defined as the remainder after subtracting
   * <code>a</code> multiplied by the mathematical integer that is less than
   * or equal to this <code>Real</code> divided by <code>a</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = this -
   *           a*Math.{@link Math#floor(double) floor}(this/a);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     27
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument.
   */
  public void mod(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0)) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      if ((a.exponent == 0 && a.mantissa == 0))
        makeNan();
      else
        sign = a.sign;
      return;
    }
    if ((a.exponent < 0 && a.mantissa == 0)) {
      if (sign != a.sign)
        makeInfinity(a.sign);
      return;
    }
    if ((a.exponent == 0 && a.mantissa == 0)) {
      makeZero(a.sign);
      return;
    }
    modInternal(a,0);
  }

  /**
   * Calculates the logical <i>AND</i> of this <code>Real</code> and
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>Semantics of bitwise logical operations exactly mimic those of
   * Java's bitwise integer operators. In these operations, the
   * internal binary representation of the numbers are used. If the
   * values represented by the operands are not mathematical
   * integers, the fractional bits are also included in the operation.
   *
   * <p>Negative numbers are interpreted as two's-complement,
   * generalized to real numbers: Negating the number inverts all
   * bits, including an infinite number of 1-bits before the radix
   * point and an infinite number of 1-bits after the radix point. The
   * infinite number of 1-bits after the radix is rounded upwards
   * producing an infinite number of 0-bits, until the first 0-bit is
   * encountered which will be switched to a 1 (rounded or not, these
   * two forms are mathematically equivalent). For example, the number
   * "1" negated, becomes (in binary form)
   * <code>...1111110.111111....</code> Rounding of the infinite
   * number of 1's after the radix gives the number
   * <code>...1111111.000000...</code>, which is exactly the way we
   * usually see "-1" as two's-complement.
   *
   * <p>This method calculates a negative value if and only
   * if this and <code>a</code> are both negative.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>int</code><i>&nbsp;code:</i></td><td>
   *     <code>this &= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.5
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument
   */
  public void and(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      makeZero();
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      if (!(this.exponent < 0 && this.mantissa == 0) && (this.sign!=0)) {
        { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
      } else if (!(a.exponent < 0 && a.mantissa == 0) && (a.sign!=0))
        ; // ASSIGN(this,this)
      else if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) &&
               (this.sign!=0) && (a.sign!=0))
        ; // makeInfinity(1)
      else
        makeZero();
      return;
    }

    byte s;
    int e;
    long m;
    if (exponent >= a.exponent) {
      s = a.sign;
      e = a.exponent;
      m = a.mantissa;
    } else {
      s = sign;
      e = exponent;
      m = mantissa;
      sign = a.sign;
      exponent = a.exponent;
      mantissa = a.mantissa;
    }
    int shift = exponent-e;
    if (shift>=64) {
      if (s == 0)
        makeZero(sign);
      return;
    }
    if (s != 0)
      m = -m;
    if ((this.sign!=0))
      mantissa = -mantissa;

    mantissa &= m>>shift;
    sign = 0;
    if (mantissa < 0) {
      mantissa = -mantissa;
      sign = 1;
    }
    normalize();
  }

  /**
   * Calculates the logical <i>OR</i> of this <code>Real</code> and
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>See {@link #and(Real)} for an explanation of the
   * interpretation of a <code>Real</code> in bitwise operations.
   * This method calculates a negative value if and only
   * if either this or <code>a</code> is negative.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>int</code><i>&nbsp;code:</i></td><td>
   *     <code>this |= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.6
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument
   */
  public void or(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      if ((this.exponent == 0 && this.mantissa == 0))
        { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      if (!(this.exponent < 0 && this.mantissa == 0) && (this.sign!=0))
        ; // ASSIGN(this,this);
      else if (!(a.exponent < 0 && a.mantissa == 0) && (a.sign!=0)) {
        { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
      } else
        makeInfinity(sign | a.sign);
      return;
    }

    byte s;
    int e;
    long m;
    if (((this.sign!=0) && exponent <= a.exponent) ||
        ((a.sign==0) && exponent >= a.exponent))
    {
      s = a.sign;
      e = a.exponent;
      m = a.mantissa;
    } else {
      s = sign;
      e = exponent;
      m = mantissa;
      sign = a.sign;
      exponent = a.exponent;
      mantissa = a.mantissa;
    }
    int shift = exponent-e;
    if (shift>=64 || shift<=-64)
      return;

    if (s != 0)
      m = -m;
    if ((this.sign!=0))
      mantissa = -mantissa;

    if (shift>=0)
      mantissa |= m>>shift;
    else
      mantissa |= m<<(-shift);

    sign = 0;
    if (mantissa < 0) {
      mantissa = -mantissa;
      sign = 1;
    }
    normalize();
  }

  /**
   * Calculates the logical <i>XOR</i> of this <code>Real</code> and
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>See {@link #and(Real)} for an explanation of the
   * interpretation of a <code>Real</code> in bitwise operations.
   * This method calculates a negative value if and only
   * if exactly one of this and <code>a</code> is negative.
   *
   * <p>The operation <i>NOT</i> has been omitted in this library
   * because it cannot be generalized to fractional numbers. If this
   * <code>Real</code> represents a mathematical integer, the
   * operation <i>NOT</i> can be calculated as "this <i>XOR</i> -1",
   * which is equivalent to "this <i>XOR</i>
   * <code>/FFFFFFFF.0000</code>".
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>int</code><i>&nbsp;code:</i></td><td>
   *     <code>this ^= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.5
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument
   */
  public void xor(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
      if ((this.exponent == 0 && this.mantissa == 0))
        { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      makeInfinity(sign ^ a.sign);
      return;
    }

    byte s;
    int e;
    long m;
    if (exponent >= a.exponent) {
      s = a.sign;
      e = a.exponent;
      m = a.mantissa;
    } else {
      s = sign;
      e = exponent;
      m = mantissa;
      sign = a.sign;
      exponent = a.exponent;
      mantissa = a.mantissa;
    }
    int shift = exponent-e;
    if (shift>=64)
      return;

    if (s != 0)
      m = -m;
    if ((this.sign!=0))
      mantissa = -mantissa;

    mantissa ^= m>>shift;

    sign = 0;
    if (mantissa < 0) {
      mantissa = -mantissa;
      sign = 1;
    }
    normalize();
  }

  /**
   * Calculates the value of this <code>Real</code> <i>AND NOT</i>
   * <code>a</code>. The opeation is read as "bit clear".
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>See {@link #and(Real)} for an explanation of the
   * interpretation of a <code>Real</code> in bitwise operations.
   * This method calculates a negative value if and only
   * if this is negative and not <code>a</code> is negative.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>int</code><i>&nbsp;code:</i></td><td>
   *     <code>this &= ~a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.5
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument
   */
  public void bic(Real a) {
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0))
      return;
    if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
      if (!(this.exponent < 0 && this.mantissa == 0)) {
        if ((this.sign!=0))
          if ((a.sign!=0))
            makeInfinity(0);
          else
            makeInfinity(1);
      } else if ((a.sign!=0)) {
        if ((a.exponent < 0 && a.mantissa == 0))
          makeInfinity(0);
        else
          makeZero();
      }
      return;
    }

    int shift = exponent-a.exponent;
    if (shift>=64 || (shift<=-64 && (this.sign==0)))
      return;

    long m = a.mantissa;
    if ((a.sign!=0))
      m = -m;
    if ((this.sign!=0))
      mantissa = -mantissa;

    if (shift<0) {
      if ((this.sign!=0)) {
        if (shift<=-64)
          mantissa = ~m;
        else
          mantissa = (mantissa>>(-shift)) & ~m;
        exponent = a.exponent;
      } else
        mantissa &= ~(m<<(-shift));
    } else
      mantissa &= ~(m>>shift);

    sign = 0;
    if (mantissa < 0) {
      mantissa = -mantissa;
      sign = 1;
    }
    normalize();
  }

  // Overloaded metohds, int arguments

  /**
   * Calculates the sum of this <code>Real</code> and the integer
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this += a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.8
   * </td></tr></table>
   *
   * @param a the <code>int</code> to add to this.
   */
  public Real add(int a) {
    tmp0.assign(a);
    add(tmp0);
    return this;
  }

  /**
   * Calculates the difference between this <code>Real</code> and the
   * integer <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this -= a;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     2.4
   * </td></tr></table>
   *
   * @param a the <code>int</code> to subtract from this.
   */
  public void sub(int a) {
    tmp0.assign(a);
    sub(tmp0);
  }

  /**
   * Calculates the quotient of the integer <code>a</code> and this
   * <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = a/this;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     3.9
   * </td></tr></table>
   *
   * @param a the <code>int</code> to be divided by this.
   */
  public void rdiv(int a) {
    tmp0.assign(a);
    rdiv(tmp0);
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is equal to
   * the integer <code>a</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this == a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.7
   * </td></tr></table>
   *
   * @param a the <code>int</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is
   *     equal to the integer <code>a</code>. <code>false</code>
   *     otherwise.
   */
  public boolean equalTo(int a) {
    tmp0.assign(a);
    return equalTo(tmp0);
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is not equal to
   * the integer <code>a</code>.
   * If this <code>Real</code> is NaN, <code>false</code> is always
   * returned.
   * This distinguishes <code>notEqualTo(a)</code> from the expression
   * <code>!equalTo(a)</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this != a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.7
   * </td></tr></table>
   *
   * @param a the <code>int</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is not
   *     equal to the integer <code>a</code>. <code>false</code>
   *     otherwise, or if this <code>Real</code> is NaN.
   */
  public boolean notEqualTo(int a) {
    tmp0.assign(a);
    return notEqualTo(tmp0);
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is less than
   * the integer <code>a</code>.
   * If this <code>Real</code> is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &lt; a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.7
   * </td></tr></table>
   *
   * @param a the <code>int</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is less
   *     than the integer <code>a</code>. <code>false</code>
   *     otherwise, or if this <code>Real</code> is NaN.
   */
  public boolean lessThan(int a) {
    tmp0.assign(a);
    return lessThan(tmp0);
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is less than or
   * equal to the integer <code>a</code>.
   * If this <code>Real</code> is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &lt;= a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.7
   * </td></tr></table>
   *
   * @param a the <code>int</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is less
   *     than or equal to the integer <code>a</code>.
   *     <code>false</code> otherwise, or if this <code>Real</code> is NaN.
   */
  public boolean lessEqual(int a) {
    tmp0.assign(a);
    return lessEqual(tmp0);
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is greater than
   * the integer <code>a</code>.
   * If this <code>Real</code> is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &gt; a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.7
   * </td></tr></table>
   *
   * @param a the <code>int</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is
   *     greater than the integer <code>a</code>.
   *     <code>false</code> otherwise, or if this <code>Real</code> is NaN.
   */
  public boolean greaterThan(int a) {
    tmp0.assign(a);
    return greaterThan(tmp0);
  }

  /**
   * Returns <code>true</code> if this <code>Real</code> is greater than
   * or equal to the integer <code>a</code>.
   * If this <code>Real</code> is NaN, <code>false</code> is always
   * returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>(this &gt;= a)</code>
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     1.7
   * </td></tr></table>
   *
   * @param a the <code>int</code> to compare to this.
   * @return <code>true</code> if the value represented by this object is
   *     greater than or equal to the integer <code>a</code>.
   *     <code>false</code> otherwise, or if this <code>Real</code> is NaN.
   */
  public boolean greaterEqual(int a) {
    tmp0.assign(a);
    return greaterEqual(tmp0);
  }

  private int compare(int a) {
    tmp0.assign(a);
    return compare(tmp0);
  }

  /**
   * Calculates the square root of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#sqrt(double) sqrt}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     19
   * </td></tr></table>
   */
  public void sqrt() {
    /*
     * Adapted from:
     * Cephes Math Library Release 2.2:  December, 1990
     * Copyright 1984, 1990 by Stephen L. Moshier
     *
     * sqrtl.c
     *
     * long double sqrtl(long double x);
     */
    if ((this.exponent < 0 && this.mantissa != 0))
      return;
    if ((this.exponent == 0 && this.mantissa == 0)) {
      sign=0;
      return;
    }
    if ((this.sign!=0)) {
      makeNan();
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0))
      return;

    // Save X
    { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; };

    // normalize to range [0.5, 1)
    int e = exponent-0x3fffffff;
    exponent = 0x3fffffff;

    // quadratic approximation, relative error 6.45e-4
    { recipTmp2.mantissa = this.mantissa; recipTmp2.exponent = this.exponent; recipTmp2.sign = this.sign; };
    { sqrtTmp.sign=(byte)1; sqrtTmp.exponent=0x3ffffffd; sqrtTmp.mantissa=0x68a7e193370ff21bL; };//-0.2044058315473477195990
    mul(sqrtTmp);
    { sqrtTmp.sign=(byte)0; sqrtTmp.exponent=0x3fffffff; sqrtTmp.mantissa=0x71f1e120690deae8L; };//0.89019407351052789754347
    add(sqrtTmp);
    mul(recipTmp2);
    { sqrtTmp.sign=(byte)0; sqrtTmp.exponent=0x3ffffffe; sqrtTmp.mantissa=0x5045ee6baf28677aL; };//0.31356706742295303132394
    add(sqrtTmp);

    // adjust for odd powers of 2
    if ((e&1) != 0)
      mul(SQRT2);

    // calculate exponent
    exponent += e>>1;

    // Newton iteratios:
    //   Yn+1 = (Yn + X/Yn)/2
    for (int i=0; i<3; i++) {
      { recipTmp2.mantissa = recipTmp.mantissa; recipTmp2.exponent = recipTmp.exponent; recipTmp2.sign = recipTmp.sign; };
      recipTmp2.div(this);
      add(recipTmp2);
      scalbn(-1);
    }
  }

  /**
   * Calculates the reciprocal square root of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = 1/Math.{@link Math#sqrt(double) sqrt}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     21
   * </td></tr></table>
   */
  public void rsqrt() {
    sqrt();
    recip();
  }

  /**
   * Calculates the cube root of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The cube root of a negative value is the negative of the cube
   * root of that value's magnitude.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#cbrt(double) cbrt}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     32
   * </td></tr></table>
   */
  public void cbrt() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    byte s = sign;
    sign = 0;
    // Calculates recipocal cube root of normalized Real,
    // not zero, nan or infinity
    final long start = 0x5120000000000000L;

    // Save -A
    { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; };
    recipTmp.neg();

    // First establish approximate result

    mantissa = start-(mantissa>>>2);
    int expRmd = exponent==0 ? 2 : (exponent-1)%3;
    exponent = 0x40000000-(exponent-0x40000000-expRmd)/3;
    normalize();
    if (expRmd>0) {
      { recipTmp2.sign=(byte)0; recipTmp2.exponent=0x3fffffff; recipTmp2.mantissa=0x6597fa94f5b8f20bL; }; // cbrt(1/2)
      mul(recipTmp2);
      if (expRmd>1)
        mul(recipTmp2);
    }

    // Now perform Newton-Raphson iteration
    // Xn+1 = (4*Xn - A*Xn**4)/3

    for (int i=0; i<4; i++) {
      { recipTmp2.mantissa = this.mantissa; recipTmp2.exponent = this.exponent; recipTmp2.sign = this.sign; };
      sqr();
      sqr();
      mul(recipTmp);
      recipTmp2.scalbn(2);
      add(recipTmp2);
      mul(THIRD);
    }
    recip();
    if (!(this.exponent < 0 && this.mantissa != 0))
      sign = s;
  }

  /**
   * Calculates the n'th root of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * For odd integer n, the n'th root of a negative value is the
   * negative of the n'th root of that value's magnitude.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#pow(double,double)
   *           pow}(this,1/a);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     110
   * </td></tr></table>
   *
   * @param n the <code>Real</code> argument.
   */
  public void nroot(Real n) {
    if ((n.exponent < 0 && n.mantissa != 0)) {
      makeNan();
      return;
    }
    if (n.compare(THREE)==0) {
      cbrt(); // Most probable application of nroot...
      return;
    } else if (n.compare(TWO)==0) {
      sqrt(); // Also possible, should be optimized like this
      return;
    }
    boolean negative = false;
    if ((this.sign!=0) && n.isIntegral() && n.isOdd()) {
      negative = true;
      abs();
    }
    { tmp2.mantissa = n.mantissa; tmp2.exponent = n.exponent; tmp2.sign = n.sign; }; // Copy to temporary location in case of x.nroot(x)
    tmp2.recip();
    pow(tmp2);
    if (negative)
      neg();
  }

  /**
   * Calculates <code>sqrt(this*this+a*a)</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#hypot(double,double)
   *           hypot}(this,a);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     24
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument.
   */
  public void hypot(Real a) {
    { tmp1.mantissa = a.mantissa; tmp1.exponent = a.exponent; tmp1.sign = a.sign; }; // Copy to temporary location in case of x.hypot(x)
    tmp1.sqr();
    sqr();
    add(tmp1);
    sqrt();
  }

  private void exp2Internal(long extra) {
    if ((this.exponent < 0 && this.mantissa != 0))
      return;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      if ((this.sign!=0))
        makeZero(0);
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
      return;
    }

    // Extract integer part
    { expTmp.mantissa = this.mantissa; expTmp.exponent = this.exponent; expTmp.sign = this.sign; };
    expTmp.add(HALF);
    expTmp.floor();
    int exp = expTmp.toInteger();
    if (exp > 0x40000000) {
      makeInfinity(sign);
      return;
    }
    if (exp < -0x40000000) {
      makeZero(sign);
      return;
    }

    // Subtract integer part (this is where we need the extra accuracy)
    expTmp.neg();
    add128(extra,expTmp,0);

    /*
     * Adapted from:
     * Cephes Math Library Release 2.7:  May, 1998
     * Copyright 1984, 1991, 1998 by Stephen L. Moshier
     *
     * exp2l.c
     *
     * long double exp2l(long double x);
     */

    // Now -0.5<X<0.5
    // rational approximation
    // exp2(x) = 1 + 2x P(x�)/(Q(x�) - x P(x�))
    { expTmp2.mantissa = this.mantissa; expTmp2.exponent = this.exponent; expTmp2.sign = this.sign; };
    expTmp2.sqr();

    // P(x�)
    { expTmp.sign=(byte)0; expTmp.exponent=0x40000005; expTmp.mantissa=0x793ace15b56b7fecL; };//60.614853552242266094567
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x4000000e; expTmp3.mantissa=0x764ef8cf96e29a13L; };//30286.971917562792508623
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000014; expTmp3.mantissa=0x7efa0173e820bf60L; };//2080384.3631901852422887
    expTmp.add(expTmp3);

    mul(expTmp);

    // Q(x�)
    expTmp.assign(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x4000000a; expTmp3.mantissa=0x6d549a6b4dc9abadL; };//1749.2876999891839021063
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000012; expTmp3.mantissa=0x5002d27836ba71c6L; };//327725.15434906797273099
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000016; expTmp3.mantissa=0x5b98206867dd59bfL; };//6002720.4078348487957118
    expTmp.add(expTmp3);

    expTmp.sub(this);
    div(expTmp);
    scalbn(1);
    add(ONE);

    // Scale by power of 2
    scalbn(exp);
  }

  /**
   * Calculates <i>e</i> raised to the power of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#exp(double) exp}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     31
   * </td></tr></table>
   */
  public void exp() {
    { expTmp.sign=(byte)0; expTmp.exponent=0x40000000; expTmp.mantissa=0x5c551d94ae0bf85dL; }; // log2(e)
    long extra = mul128(0,expTmp,0xdf43ff68348e9f44L);
    exp2Internal(extra);
  }

  /**
   * Calculates 2 raised to the power of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#exp(double) exp}(this *
   *           Math.{@link Math#log(double) log}(2));</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     27
   * </td></tr></table>
   */
  public void exp2() {
    exp2Internal(0);
  }

  /**
   * Calculates 10 raised to the power of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#exp(double) exp}(this *
   *           Math.{@link Math#log(double) log}(10));</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     31
   * </td></tr></table>
   */
  public void exp10() {
    { expTmp.sign=(byte)0; expTmp.exponent=0x40000001; expTmp.mantissa=0x6a4d3c25e68dc57fL; }; // log2(10)
    long extra = mul128(0,expTmp,0x2495fb7fa6d7eda6L);
    exp2Internal(extra);
  }

  private int lnInternal()
  {
    if ((this.exponent < 0 && this.mantissa != 0))
      return 0;
    if ((this.sign!=0)) {
      makeNan();
      return 0;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      makeInfinity(1);
      return 0;
    }
    if ((this.exponent < 0 && this.mantissa == 0))
      return 0;

    /*
     * Adapted from:
     * Cephes Math Library Release 2.7:  May, 1998
     * Copyright 1984, 1990, 1998 by Stephen L. Moshier
     *
     * logl.c
     *
     * long double logl(long double x);
     */

    // normalize to range [0.5, 1)
    int e = exponent-0x3fffffff;
    exponent = 0x3fffffff;

    // rational appriximation
    // log(1+x) = x - x�/2 + x� P(x)/Q(x)
    if (this.compare(SQRT1_2) < 0) {
      e--;
      exponent++;
    }
    sub(ONE);
    { expTmp2.mantissa = this.mantissa; expTmp2.exponent = this.exponent; expTmp2.sign = this.sign; };

    // P(x)
    { this.sign=(byte)0; this.exponent=0x3ffffff1; this.mantissa=0x5ef0258ace5728ddL; };//4.5270000862445199635215E-5
    mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x3ffffffe; expTmp3.mantissa=0x7fa06283f86a0ce8L; };//0.4985410282319337597221
    add(expTmp3);
    mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000002; expTmp3.mantissa=0x69427d1bd3e94ca1L; };//6.5787325942061044846969
    add(expTmp3);
    mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000004; expTmp3.mantissa=0x77a5ce2e32e7256eL; };//29.911919328553073277375
    add(expTmp3);
    mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000005; expTmp3.mantissa=0x79e63ae1b0cd4222L; };//60.949667980987787057556
    add(expTmp3);
    mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000005; expTmp3.mantissa=0x7239d65d1e6840d6L; };//57.112963590585538103336
    add(expTmp3);
    mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000004; expTmp3.mantissa=0x502880b6660c265fL; };//20.039553499201281259648
    add(expTmp3);

    // Q(x)
    { expTmp.mantissa = expTmp2.mantissa; expTmp.exponent = expTmp2.exponent; expTmp.sign = expTmp2.sign; };
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000003; expTmp3.mantissa=0x7880d67a40f8dc5cL; };//15.062909083469192043167
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000006; expTmp3.mantissa=0x530c2d4884d25e18L; };//83.047565967967209469434
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000007; expTmp3.mantissa=0x6ee19643f3ed5776L; };//221.76239823732856465394
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000008; expTmp3.mantissa=0x4d465177242295efL; };//309.09872225312059774938
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000007; expTmp3.mantissa=0x6c36c4f923819890L; };//216.42788614495947685003
    expTmp.add(expTmp3);
    expTmp.mul(expTmp2);
    { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000005; expTmp3.mantissa=0x783cc111991239a3L; };//60.118660497603843919306
    expTmp.add(expTmp3);

    div(expTmp);
    { expTmp3.mantissa = expTmp2.mantissa; expTmp3.exponent = expTmp2.exponent; expTmp3.sign = expTmp2.sign; };
    expTmp3.sqr();
    mul(expTmp3);
    mul(expTmp2);
    expTmp3.scalbn(-1);
    sub(expTmp3);
    add(expTmp2);

    return e;
  }

  /**
   * Calculates the natural logarithm (base-<i>e</i>) of this
   * <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#log(double) log}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     51
   * </td></tr></table>
   */
  public void ln() {
    int exp = lnInternal();

    expTmp.assign(exp);
    expTmp.mul(LN2);
    add(expTmp);
  }

  /**
   * Calculates the base-2 logarithm of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#log(double) log}(this)/Math.{@link
   *           Math#log(double) log}(2);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     51
   * </td></tr></table>
   */
  public void log2() {
    int exp = lnInternal();

    mul(LOG2E);
    add(exp);
  }

  /**
   * Calculates the base-10 logarithm of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#log10(double) log10}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     53
   * </td></tr></table>
   */
  public void log10() {
    int exp = lnInternal();

    expTmp.assign(exp);
    expTmp.mul(LN2);
    add(expTmp);
    mul(LOG10E);
  }

  /**
   * Calculates the closest power of 10 that is less than or equal to this
   * <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The base-10 exponent of the result is returned.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>int exp = (int)(Math.{@link Math#floor(double)
   *       floor}(Math.{@link Math#log10(double) log10}(this)));
   *       <br>this = Math.{@link Math#pow(double,double) pow}(10, exp);<br>
   *           return exp;</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     3.6
   * </td></tr></table>
   *
   * @return the base-10 exponent
   */
  public int lowPow10() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return 0;
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    // Approximate log10 using exponent only
    int e = exponent - 0x40000000;
    if (e<0) // it's important to achieve floor(exponent*ln2/ln10)
      e = -(int)(((-e)*0x4d104d43L+((1L<<32)-1)) >> 32);
    else
      e = (int)(e*0x4d104d43L >> 32);
    // Now, e < log10(this) < e+1
    { this.mantissa = TEN.mantissa; this.exponent = TEN.exponent; this.sign = TEN.sign; };
    pow(e);
    { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; };
    tmp3.mul10();
    if (tmp3.compare(tmp2) <= 0) {
      // First estimate of log10 was too low
      e++;
      { this.mantissa = tmp3.mantissa; this.exponent = tmp3.exponent; this.sign = tmp3.sign; };
    }
    return e;
  }

  /**
   * Calculates the value of this <code>Real</code> raised to the power of
   * the integer <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#pow(double,double) pow}(this, a);</code>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     84
   * </td></tr></table>
   */
  public void pow(int a) {
    // Calculate power of integer by successive squaring
    boolean recp=false;
    if (a < 0) {
      a = -a; // Also works for 0x80000000
      recp = true;
    }
    long extra = 0, expTmpExtra = 0;
    { expTmp.mantissa = this.mantissa; expTmp.exponent = this.exponent; expTmp.sign = this.sign; };
    { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
    for (; a!=0; a>>>=1) {
      if ((a & 1) != 0)
        extra = mul128(extra,expTmp,expTmpExtra);
      expTmpExtra = expTmp.mul128(expTmpExtra,expTmp,expTmpExtra);
    }
    if (recp)
      extra = recip128(extra);
    roundFrom128(extra);
  }

  /**
   * Calculates the value of this <code>Real</code> raised to the power of
   * <code>a</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p> Special cases:
   * <ul>
   * <li> if a is 0.0 or -0.0 then result is 1.0
   * <li> if a is NaN then result is NaN
   * <li> if this is NaN and a is not zero then result is NaN
   * <li> if a is 1.0 then result is this
   * <li> if |this| > 1.0 and a is +Infinity then result is +Infinity
   * <li> if |this| < 1.0 and a is -Infinity then result is +Infinity
   * <li> if |this| > 1.0 and a is -Infinity then result is +0
   * <li> if |this| < 1.0 and a is +Infinity then result is +0
   * <li> if |this| = 1.0 and a is �Infinity then result is NaN
   * <li> if this = +0 and a > 0 then result is +0
   * <li> if this = +0 and a < 0 then result is +Inf
   * <li> if this = -0 and a > 0, and odd integer then result is -0
   * <li> if this = -0 and a < 0, and odd integer then result is -Inf
   * <li> if this = -0 and a > 0, not odd integer then result is +0
   * <li> if this = -0 and a < 0, not odd integer then result is +Inf
   * <li> if this = +Inf and a > 0 then result is +Inf
   * <li> if this = +Inf and a < 0 then result is +0
   * <li> if this = -Inf and a not integer then result is NaN
   * <li> if this = -Inf and a > 0, and odd integer then result is -Inf
   * <li> if this = -Inf and a > 0, not odd integer then result is +Inf
   * <li> if this = -Inf and a < 0, and odd integer then result is -0
   * <li> if this = -Inf and a < 0, not odd integer then result is +0
   * <li> if this < 0 and a not integer then result is NaN
   * <li> if this < 0 and a odd integer then result is -(|this|<sup>a</sup>)
   * <li> if this < 0 and a not odd integer then result is |this|<sup>a</sup>
   * <li> else result is exp(ln(this)*a)
   * </ul>
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#pow(double,double) pow}(this, a);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     110
   * </td></tr></table>
   *
   * @param a the <code>Real</code> argument.
   */
  public void pow(Real a) {
    if ((a.exponent == 0 && a.mantissa == 0)) {
      { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
      return;
    }
    if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
      makeNan();
      return;
    }
    if (a.compare(ONE)==0)
      return;
    if ((a.exponent < 0 && a.mantissa == 0)) {
      { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
      tmp1.abs();
      int test = tmp1.compare(ONE);
      if (test>0) {
        if ((a.sign==0))
          makeInfinity(0);
        else
          makeZero();
      } else if (test<0) {
        if ((a.sign!=0))
          makeInfinity(0);
        else
          makeZero();
      } else {
        makeNan();
      }
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      if ((this.sign==0)) {
        if ((a.sign==0))
          makeZero();
        else
          makeInfinity(0);
      } else {
        if (a.isIntegral() && a.isOdd()) {
          if ((a.sign==0))
            makeZero(1);
          else
            makeInfinity(1);
        } else {
          if ((a.sign==0))
            makeZero();
          else
            makeInfinity(0);
        }
      }
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0)) {
      if ((this.sign==0)) {
        if ((a.sign==0))
          makeInfinity(0);
        else
          makeZero();
      } else {
        if (a.isIntegral()) {
          if (a.isOdd()) {
            if ((a.sign==0))
              makeInfinity(1);
            else
              makeZero(1);
          } else {
            if ((a.sign==0))
              makeInfinity(0);
            else
              makeZero();
          }
        } else {
          makeNan();
        }
      }
      return;
    }

    if (a.isIntegral() && a.exponent <= 0x4000001e) {
      pow(a.toInteger());
      return;
    }

    byte s=0;
    if ((this.sign!=0)) {
      if (a.isIntegral()) {
        if (a.isOdd())
          s = 1;
      } else {
        makeNan();
        return;
      }
      sign = 0;
    }

    { tmp1.mantissa = a.mantissa; tmp1.exponent = a.exponent; tmp1.sign = a.sign; };
    if (tmp1.exponent <= 0x4000001e) {
      // For increased accuracy, exponentiate with integer part of
      // exponent by successive squaring
      // (I really don't know why this works)
      { tmp2.mantissa = tmp1.mantissa; tmp2.exponent = tmp1.exponent; tmp2.sign = tmp1.sign; };
      tmp2.floor();
      { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; };
      tmp3.pow(tmp2.toInteger());
      tmp1.sub(tmp2);
    } else {
      { tmp3.mantissa = ONE.mantissa; tmp3.exponent = ONE.exponent; tmp3.sign = ONE.sign; };
    }

    // Do log2 and maintain accuracy
    int e = lnInternal();
    { tmp2.sign=(byte)0; tmp2.exponent=0x40000000; tmp2.mantissa=0x5c551d94ae0bf85dL; }; // log2(e)
    long extra = mul128(0,tmp2,0xdf43ff68348e9f44L);
    tmp2.assign(e);
    extra = add128(extra,tmp2,0);

    // Do exp2 of this multiplied by (fractional part of) exponent
    extra = tmp1.mul128(0,this,extra);
    tmp1.exp2Internal(extra);
    { this.mantissa = tmp1.mantissa; this.exponent = tmp1.exponent; this.sign = tmp1.sign; };
    mul(tmp3);

    if (!(this.exponent < 0 && this.mantissa != 0))
      sign = s;
  }

  private void sinInternal() {
    /*
     * Adapted from:
     * Cephes Math Library Release 2.7:  May, 1998
     * Copyright 1985, 1990, 1998 by Stephen L. Moshier
     *
     * sinl.c
     *
     * long double sinl(long double x);
     */

    // X<PI/4
    // polynomial approximation
    // sin(x) = x + x� P(x�)
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    tmp2.sqr();
    { this.sign=(byte)1; this.exponent=0x3fffffd7; this.mantissa=0x6aa891c4f0eb2713L; };//-7.578540409484280575629E-13
    mul(tmp2);
    { tmp3.sign=(byte)0; tmp3.exponent=0x3fffffdf; tmp3.mantissa=0x58482311f383326cL; };//1.6058363167320443249231E-10
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3fffffe6; tmp3.mantissa=0x6b9914a35f9a00d8L; };//-2.5052104881870868784055E-8
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)0; tmp3.exponent=0x3fffffed; tmp3.mantissa=0x5c778e94cc22e47bL; };//2.7557319214064922217861E-6
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffff3; tmp3.mantissa=0x680680680629b28aL; };//-1.9841269841254799668344E-4
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)0; tmp3.exponent=0x3ffffff9; tmp3.mantissa=0x4444444444442b4dL; };//8.3333333333333225058715E-3
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffd; tmp3.mantissa=0x555555555555554cL; };//-1.6666666666666666640255E-1
    add(tmp3);
    mul(tmp2);
    mul(tmp1);
    add(tmp1);
  }

  private void cosInternal() {
    /*
     * Adapted from:
     * Cephes Math Library Release 2.7:  May, 1998
     * Copyright 1985, 1990, 1998 by Stephen L. Moshier
     *
     * sinl.c
     *
     * long double cosl(long double x);
     */

    // X<PI/4
    // polynomial approximation
    // cos(x) = 1 - x�/2 + x**4 Q(x�)
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    tmp2.sqr();
    { this.sign=(byte)0; this.exponent=0x3fffffd3; this.mantissa=0x6aaf461d37ccba1bL; };//4.7377507964246204691685E-14
    mul(tmp2);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3fffffdb; tmp3.mantissa=0x64e4c907ac7a179bL; };//-1.147028484342535976567E-11
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)0; tmp3.exponent=0x3fffffe3; tmp3.mantissa=0x47bb632432cf29a8L; };//2.0876754287081521758361E-9
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3fffffea; tmp3.mantissa=0x49f93edd7ae32696L; };//-2.7557319214999787979814E-7
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)0; tmp3.exponent=0x3ffffff0; tmp3.mantissa=0x68068068063329f7L; };//2.4801587301570552304991E-5L
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffff6; tmp3.mantissa=0x5b05b05b05b03db3L; };//-1.3888888888888872993737E-3
    add(tmp3);
    mul(tmp2);
    { tmp3.sign=(byte)0; tmp3.exponent=0x3ffffffb; tmp3.mantissa=0x555555555555554dL; };//4.1666666666666666609054E-2
    add(tmp3);
    mul(tmp2);
    sub(HALF);
    mul(tmp2);
    add(ONE);
  }

  /**
   * Calculates the trigonometric sine of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The input value is treated as an angle measured in radians.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#sin(double) sin}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     28
   * </td></tr></table>
   */
  public void sin() {
    if (!(this.exponent >= 0 && this.mantissa != 0)) {
      if (!(this.exponent == 0 && this.mantissa == 0))
        makeNan();
      return;
    }

    // Since sin(-x) = -sin(x) we can make sure that x > 0
    boolean negative = false;
    if ((this.sign!=0)) {
      abs();
      negative = true;
    }

    // Then reduce the argument to the range of 0 < x < pi*2
    if (this.compare(PI2) > 0)
      modInternal(PI2,0x62633145c06e0e69L);

    // Since sin(pi*2 - x) = -sin(x) we can reduce the range 0 < x < pi
    if (this.compare(PI) > 0) {
      sub(PI2);
      neg();
      negative = !negative;
    }

    // Since sin(x) = sin(pi - x) we can reduce the range to 0 < x < pi/2
    if (this.compare(PI_2) > 0) {
      sub(PI);
      neg();
    }

    // Since sin(x) = cos(pi/2 - x) we can reduce the range to 0 < x < pi/4
    if (this.compare(PI_4) > 0) {
      sub(PI_2);
      neg();
      cosInternal();
    } else {
      sinInternal();
    }

    if (negative)
      neg();
    if ((this.exponent == 0 && this.mantissa == 0))
      abs(); // Remove confusing "-"
  }

  /**
   * Calculates the trigonometric cosine of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The input value is treated as an angle measured in radians.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#cos(double) cos}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     1 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     37
   * </td></tr></table>
   */
  public void cos() {
    if ((this.exponent == 0 && this.mantissa == 0)) {
      { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
      return;
    }
    if ((this.sign!=0))
      abs();
    if (this.compare(PI_4) < 0) {
      cosInternal();
    } else {
      add(PI_2);
      sin();
    }
  }

  /**
   * Calculates the trigonometric tangent of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The input value is treated as an angle measured in radians.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#tan(double) tan}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     70
   * </td></tr></table>
   */
  public void tan() {
    { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; };
    tmp4.cos();
    sin();
    div(tmp4);
  }

  /**
   * Calculates the trigonometric arc sine of this <code>Real</code>,
   * in the range -&pi;/2 to &pi;/2.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#asin(double) asin}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     3 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     68
   * </td></tr></table>
   */
  public void asin() {
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    sqr();
    neg();
    add(ONE);
    rsqrt();
    mul(tmp1);
    atan();
  }

  /**
   * Calculates the trigonometric arc cosine of this <code>Real</code>,
   * in the range 0.0 to &pi;.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#acos(double) acos}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     67
   * </td></tr></table>
   */
  public void acos() {
    boolean negative = (this.sign!=0);
    abs();

    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    sqr();
    neg();
    add(ONE);
    sqrt();
    div(tmp1);
    atan();

    if (negative) {
      neg();
      add(PI);
    }
  }

  /**
   * Calculates the trigonometric arc tangent of this <code>Real</code>,
   * in the range -&pi;/2 to &pi;/2.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#atan(double) atan}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     37
   * </td></tr></table>
   */
  public void atan() {
    /*
     * Adapted from:
     * Cephes Math Library Release 2.7:  May, 1998
     * Copyright 1984, 1990, 1998 by Stephen L. Moshier
     *
     * atanl.c
     *
     * long double atanl(long double x);
     */
    if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa != 0))
      return;
    if ((this.exponent < 0 && this.mantissa == 0)) {
      byte s = sign;
      { this.mantissa = PI_2.mantissa; this.exponent = PI_2.exponent; this.sign = PI_2.sign; };
      sign = s;
      return;
    }

    byte s = sign;
    sign = 0;

    // range reduction
    boolean addPI_2 = false;
    boolean addPI_4 = false;
    { tmp1.mantissa = SQRT2.mantissa; tmp1.exponent = SQRT2.exponent; tmp1.sign = SQRT2.sign; };
    tmp1.add(ONE);
    if (this.compare(tmp1) > 0) {
      addPI_2 = true;
      recip();
      neg();
    } else {
      tmp1.sub(TWO);
      if (this.compare(tmp1) > 0) {
        addPI_4 = true;
        { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
        tmp1.add(ONE);
        sub(ONE);
        div(tmp1);
      }
    }

    // Now |X|<sqrt(2)-1
    // rational approximation
    // atan(x) = x + x� P(x�)/Q(x�)
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    tmp2.sqr();
    mul(tmp2);

    { tmp3.sign=(byte)1; tmp3.exponent=0x3fffffff; tmp3.mantissa=0x6f2f89336729c767L; };//-0.8686381817809218753544
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)1; tmp4.exponent=0x40000003; tmp4.mantissa=0x7577d35fd03083f3L; };//-14.683508633175792446076
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)1; tmp4.exponent=0x40000005; tmp4.mantissa=0x7ff42abff948a9f7L; };//-63.976888655834347413154
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)1; tmp4.exponent=0x40000006; tmp4.mantissa=0x63fd1f9f76d37cebL; };//-99.988763777265819915721
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)1; tmp4.exponent=0x40000005; tmp4.mantissa=0x65c9c9b0b55e5b62L; };//-50.894116899623603312185
    tmp3.add(tmp4);

    mul(tmp3);

    { tmp3.mantissa = tmp2.mantissa; tmp3.exponent = tmp2.exponent; tmp3.sign = tmp2.sign; };
    { tmp4.sign=(byte)0; tmp4.exponent=0x40000004; tmp4.mantissa=0x5bed73b744a72a6aL; };//22.981886733594175366172
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)0; tmp4.exponent=0x40000007; tmp4.mantissa=0x47fed7d13d233b5cL; };//143.99096122250781605352
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)0; tmp4.exponent=0x40000008; tmp4.mantissa=0x5a5c35f774e071d5L; };//361.44079386152023162701
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)0; tmp4.exponent=0x40000008; tmp4.mantissa=0x61e4d84c2853d5e0L; };//391.57570175111990631099
    tmp3.add(tmp4);
    tmp3.mul(tmp2);
    { tmp4.sign=(byte)0; tmp4.exponent=0x40000007; tmp4.mantissa=0x4c5757448806c48eL; };//152.68235069887081006606
    tmp3.add(tmp4);

    div(tmp3);
    add(tmp1);

    if (addPI_2)
      add(PI_2);
    if (addPI_4)
      add(PI_4);
    if (s != 0)
      neg();
  }

  /**
   * Calculates the trigonometric arc tangent of this
   * <code>Real</code> divided by <code>x</code>, in the range -&pi;
   * to &pi;. The signs of both arguments are used to determine the
   * quadrant of the result. Replaces the contents of this
   * <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#atan2(double,double)
   *           atan2}(this,x);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     48
   * </td></tr></table>
   *
   * @param x the <code>Real</code> argument.
   */
  public void atan2(Real x) {
    if ((this.exponent < 0 && this.mantissa != 0) || (x.exponent < 0 && x.mantissa != 0) || ((this.exponent < 0 && this.mantissa == 0) && (x.exponent < 0 && x.mantissa == 0))) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0) && (x.exponent == 0 && x.mantissa == 0))
      return;

    byte s = sign;
    byte s2 = x.sign;
    sign = 0;
    x.sign = 0;

    div(x);
    atan();

    if (s2 != 0) {
      neg();
      add(PI);
    }

    sign = s;
  }

  /**
   * Calculates the hyperbolic sine of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#sinh(double) sinh}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     67
   * </td></tr></table>
   */
  public void sinh() {
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.neg();
    tmp1.exp();
    exp();
    sub(tmp1);
    scalbn(-1);
  }

  /**
   * Calculates the hyperbolic cosine of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#cosh(double) cosh}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     66
   * </td></tr></table>
   */
  public void cosh() {
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.neg();
    tmp1.exp();
    exp();
    add(tmp1);
    scalbn(-1);
  }

  /**
   * Calculates the hyperbolic tangent of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#tanh(double) tanh}(this);</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     70
   * </td></tr></table>
   */
  public void tanh() {
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.neg();
    tmp1.exp();
    exp();
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    tmp2.add(tmp1);
    sub(tmp1);
    div(tmp2);
  }

  /**
   * Calculates the hyperbolic arc sine of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     77
   * </td></tr></table>
   */
  public void asinh() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;

    // Use symmetry to prevent underflow error for very large negative values
    byte s = sign;
    sign = 0;

    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.sqr();
    tmp1.add(ONE);
    tmp1.sqrt();
    add(tmp1);
    ln();

    if (!(this.exponent < 0 && this.mantissa != 0))
      sign = s;
  }

  /**
   * Calculates the hyperbolic arc cosine of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     75
   * </td></tr></table>
   */
  public void acosh() {
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.sqr();
    tmp1.sub(ONE);
    tmp1.sqrt();
    add(tmp1);
    ln();
  }

  /**
   * Calculates the hyperbolic arc tangent of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     57
   * </td></tr></table>
   */
  public void atanh() {
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.neg();
    tmp1.add(ONE);
    add(ONE);
    div(tmp1);
    ln();
    scalbn(-1);
  }

  /**
   * Calculates the factorial of this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   * The definition is generalized to all real numbers (not only integers),
   * by using the fact that <code>(n!)={@link #gamma() gamma}(n+1)</code>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     15 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     8-190
   * </td></tr></table>
   */
  public void fact() {
    if (!(this.exponent >= 0))
      return;

    if (!this.isIntegral() || this.compare(ZERO)<0 || this.compare(200)>0)
    {
      // x<0, x>200 or not integer: fact(x) = gamma(x+1)
      add(ONE);
      gamma();
      return;
    }
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
    while (tmp1.compare(ONE) > 0) {
      mul(tmp1);
      tmp1.sub(ONE);
    }
  }

  /**
   * Calculates the gamma function for this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     100+ ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     190
   * </td></tr></table>
   */
  public void gamma() {
    if (!(this.exponent >= 0))
      return;

    // x<0: gamma(-x) = -pi/(x*gamma(x)*sin(pi*x))
    boolean negative = (this.sign!=0);
    abs();
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };

    // x<n: gamma(x) = gamma(x+m)/x*(x+1)*(x+2)*...*(x+m-1)
    // n=20
    { tmp2.mantissa = ONE.mantissa; tmp2.exponent = ONE.exponent; tmp2.sign = ONE.sign; };
    boolean divide = false;
    while (this.compare(20) < 0) {
      divide = true;
      tmp2.mul(this);
      add(ONE);
    }

    // x>n: gamma(x) = exp((x-1/2)*ln(x)-x+ln(2*pi)/2+1/12x-1/360x�+1/1260x**5
    //                     -1/1680x**7+1/1188x**9)
    { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; }; // x
    { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; };
    tmp4.sqr(); // x�

    // (x-1/2)*ln(x)-x
    ln(); { tmp5.mantissa = tmp3.mantissa; tmp5.exponent = tmp3.exponent; tmp5.sign = tmp3.sign; }; tmp5.sub(HALF); mul(tmp5); sub(tmp3);
    // + ln(2*pi)/2
    { tmp5.sign=(byte)0; tmp5.exponent=0x3fffffff; tmp5.mantissa=0x759fc72192fad29aL; }; add(tmp5);
    // + 1/12x
    tmp5.assign( 12); tmp5.mul(tmp3); tmp5.recip(); add(tmp5); tmp3.mul(tmp4);
    // - 1/360x�
    tmp5.assign( 360); tmp5.mul(tmp3); tmp5.recip(); sub(tmp5); tmp3.mul(tmp4);
    // + 1/1260x**5
    tmp5.assign(1260); tmp5.mul(tmp3); tmp5.recip(); add(tmp5); tmp3.mul(tmp4);
    // - 1/1680x**7
    tmp5.assign(1680); tmp5.mul(tmp3); tmp5.recip(); sub(tmp5); tmp3.mul(tmp4);
    // + 1/1188x**9
    tmp5.assign(1188); tmp5.mul(tmp3); tmp5.recip(); add(tmp5);

    exp();

    if (divide)
      div(tmp2);

    if (negative) {
      { tmp5.mantissa = tmp1.mantissa; tmp5.exponent = tmp1.exponent; tmp5.sign = tmp1.sign; }; // sin() uses tmp1
      // -pi/(x*gamma(x)*sin(pi*x))
      mul(tmp5);
      tmp5.scalbn(-1); tmp5.frac(); tmp5.mul(PI2); // Fixes integer inaccuracy
      tmp5.sin(); mul(tmp5); recip(); mul(PI); neg();
    }
  }

  private void erfc1Internal() {
    //                                  3       5        7        9
    //                  2     /        x       x        x        x           \
    // erfc(x) = 1 - -------- | x  -  ---  +  ----  -  ----  +  ----  -  ... |
    //               sqrt(pi) \        3      2!*5     3!*7     4!*9         /
    //
    long extra=0,tmp1Extra,tmp2Extra,tmp3Extra,tmp4Extra;
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; tmp1Extra = 0;
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    tmp2Extra = tmp2.mul128(0,tmp2,0);
    tmp2.neg();
    { tmp3.mantissa = ONE.mantissa; tmp3.exponent = ONE.exponent; tmp3.sign = ONE.sign; }; tmp3Extra = 0;
    int i=1;
    do {
      tmp1Extra = tmp1.mul128(tmp1Extra,tmp2,tmp2Extra);
      tmp4.assign(i);
      tmp3Extra = tmp3.mul128(tmp3Extra,tmp4,0);
      tmp4.assign(2*i+1);
      tmp4Extra = tmp4.mul128(0,tmp3,tmp3Extra);
      tmp4Extra = tmp4.recip128(tmp4Extra);
      tmp4Extra = tmp4.mul128(tmp4Extra,tmp1,tmp1Extra);
      extra = add128(extra,tmp4,tmp4Extra);
      i++;
    } while (exponent - tmp4.exponent < 128);
    { tmp1.sign=(byte)1; tmp1.exponent=0x40000000; tmp1.mantissa=0x48375d410a6db446L; }; // -2/sqrt(pi)
    extra = mul128(extra,tmp1,0xb8ea453fb5ff61a2L);
    extra = add128(extra,ONE,0);
    roundFrom128(extra);
  }

  private void erfc2Internal() {
    //             -x� -1
    //            e   x   /      1      3       3*5     3*5*7        \
    // erfc(x) = -------- | 1 - --- + ------ - ------ + ------ - ... |
    //           sqrt(pi) \     2x�        2        3        4       /
    //                                (2x�)    (2x�)    (2x�)

    // Calculate iteration stop criteria
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.sqr();
    { tmp2.sign=(byte)0; tmp2.exponent=0x40000000; tmp2.mantissa=0x5c3811b4bfd0c8abL; }; // 1/0.694
    tmp2.mul(tmp1);
    tmp2.sub(HALF);
    int digits = tmp2.toInteger(); // number of accurate digits = x*x/0.694-0.5
    if (digits > 64)
      digits = 64;
    tmp1.scalbn(1);
    int dxq = tmp1.toInteger()+1;

    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    recip();
    { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
    { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; };
    tmp3.sqr();
    tmp3.neg();
    tmp3.scalbn(-1);
    { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
    { tmp4.mantissa = ONE.mantissa; tmp4.exponent = ONE.exponent; tmp4.sign = ONE.sign; };
    int i=1;
    do {
      tmp4.mul(2*i-1);
      tmp4.mul(tmp3);
      add(tmp4);
      i++;
    } while (tmp4.exponent-0x40000000>-(digits+2) && 2*i-1<dxq);
    mul(tmp2);
    tmp1.sqr();
    tmp1.neg();
    tmp1.exp();
    mul(tmp1);
    { tmp1.sign=(byte)0; tmp1.exponent=0x3fffffff; tmp1.mantissa=0x48375d410a6db447L; }; // 1/sqrt(pi)
    mul(tmp1);
  }

  /**
   * Calculates the complementary error function for this <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>The complementary error function is defined as the integral from
   * x to infinity of 2/&#8730;<span style="text-decoration:
   * overline;">&pi;</span>&nbsp;�<i>e</i><sup>-t�</sup>&nbsp;dt. It is
   * related to the error function, <i>erf</i>, by the formula
   * erfc(x)=1-erf(x).
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2<sup>19</sup> ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     80-4900
   * </td></tr></table>
   */
  public void erfc() {
    if ((this.exponent < 0 && this.mantissa != 0))
      return;
    if ((this.exponent == 0 && this.mantissa == 0)) {
      { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
      return;
    }
    if ((this.exponent < 0 && this.mantissa == 0) || toInteger()>27281) {
      if ((this.sign!=0)) {
        { this.mantissa = TWO.mantissa; this.exponent = TWO.exponent; this.sign = TWO.sign; };
      } else
        makeZero(0);
      return;
    }
    byte s = sign;
    sign = 0;

    { tmp1.sign=(byte)0; tmp1.exponent=0x40000002; tmp1.mantissa=0x570a3d70a3d70a3dL; }; // 5.44
    if (this.lessThan(tmp1))
      erfc1Internal();
    else
      erfc2Internal();

    if (s != 0) {
      neg();
      add(TWO);
    }
  }

  /**
   * Calculates the inverse complementary error function for this
   * <code>Real</code>.
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     2<sup>19</sup> ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     240-5100
   * </td></tr></table>
   */
  public void inverfc() {
    if ((this.exponent < 0 && this.mantissa != 0) || (this.sign!=0) || this.greaterThan(TWO)) {
      makeNan();
      return;
    }
    if ((this.exponent == 0 && this.mantissa == 0)) {
      makeInfinity(0);
      return;
    }
    if (this.equals(TWO)) {
      makeInfinity(1);
      return;
    }
    int sign = ONE.compare(this);
    if (sign==0) {
      makeZero();
      return;
    }
    if (sign<0) {
      neg();
      add(TWO);
    }
    // Using invphi to calculate inverfc, like this
    // inverfc(x) = -invphi(x/2)/(sqrt(2))
    scalbn(-1);

    // Inverse Phi Algorithm (phi(Z)=P, so invphi(P)=Z)
    // ------------------------------------------------
    // Part 1: Numerical Approximation Method for Inverse Phi
    // This accepts input of P and outputs approximate Z as Y
    // Source:Odeh & Evans. 1974. AS 70. Applied Statistics.

    // R = sqrt(Ln(1/(Q�)))
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.ln();
    tmp1.mul(-2);
    tmp1.sqrt();

    // Y = -(R+((((P4*R+P3)*R+P2)*R+P1)*R+P0)/((((Q4*R+Q3)*R*Q2)*R+Q1)*R+Q0))
    { tmp2.sign=(byte)1; tmp2.exponent=0x3ffffff1; tmp2.mantissa=0x5f22bb0fb4698674L; }; // P4=-0.0000453642210148
    tmp2.mul(tmp1);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffa; tmp3.mantissa=0x53a731ce1ea0be15L; }; // P3=-0.0204231210245
    tmp2.add(tmp3);
    tmp2.mul(tmp1);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffe; tmp3.mantissa=0x579d2d719fc517f3L; }; // P2=-0.342242088547
    tmp2.add(tmp3);
    tmp2.mul(tmp1);
    tmp2.add(-1); // P1=-1
    tmp2.mul(tmp1);
    { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffe; tmp3.mantissa=0x527dd3193bc8dd4cL; }; // P0=-0.322232431088
    tmp2.add(tmp3);

    { tmp3.sign=(byte)0; tmp3.exponent=0x3ffffff7; tmp3.mantissa=0x7e5b0f681d161e7dL; }; // Q4=0.0038560700634
    tmp3.mul(tmp1);
    { tmp4.sign=(byte)0; tmp4.exponent=0x3ffffffc; tmp4.mantissa=0x6a05ccf9917da0a8L; }; // Q3=0.103537752850
    tmp3.add(tmp4);
    tmp3.mul(tmp1);
    { tmp4.sign=(byte)0; tmp4.exponent=0x3fffffff; tmp4.mantissa=0x43fb32c0d3c14ec4L; }; // Q2=0.531103462366
    tmp3.add(tmp4);
    tmp3.mul(tmp1);
    { tmp4.sign=(byte)0; tmp4.exponent=0x3fffffff; tmp4.mantissa=0x4b56a41226f4ba95L; }; // Q1=0.588581570495
    tmp3.add(tmp4);
    tmp3.mul(tmp1);
    { tmp4.sign=(byte)0; tmp4.exponent=0x3ffffffc; tmp4.mantissa=0x65bb9a7733dd5062L; }; // Q0=0.0993484626060
    tmp3.add(tmp4);

    tmp2.div(tmp3);
    tmp1.add(tmp2);
    tmp1.neg();
    { sqrtTmp.mantissa = tmp1.mantissa; sqrtTmp.exponent = tmp1.exponent; sqrtTmp.sign = tmp1.sign; }; // sqrtTmp and tmp5 not used by erfc() and exp()

    // Part 2: Refine to accuracy of erfc Function
    // This accepts inputs Y and P (from above) and outputs Z
    // (Using Halley's third order method for finding roots of equations)

    // Q = erfc(-Y/sqrt(2))/2-P
    { tmp5.mantissa = sqrtTmp.mantissa; tmp5.exponent = sqrtTmp.exponent; tmp5.sign = sqrtTmp.sign; };
    tmp5.mul(SQRT1_2);
    tmp5.neg();
    tmp5.erfc();
    tmp5.scalbn(-1);
    tmp5.sub(this);

    // R = Q*sqrt(2*pi)*e^(Y�/2)
    { tmp3.mantissa = sqrtTmp.mantissa; tmp3.exponent = sqrtTmp.exponent; tmp3.sign = sqrtTmp.sign; };
    tmp3.sqr();
    tmp3.scalbn(-1);
    tmp3.exp();
    tmp5.mul(tmp3);
    { tmp3.sign=(byte)0; tmp3.exponent=0x40000001; tmp3.mantissa=0x50364c7fd89c1659L; }; // sqrt(2*pi)
    tmp5.mul(tmp3);

    // Z = Y-R/(1+R*Y/2)
    { this.mantissa = sqrtTmp.mantissa; this.exponent = sqrtTmp.exponent; this.sign = sqrtTmp.sign; };
    mul(tmp5);
    scalbn(-1);
    add(ONE);
    rdiv(tmp5);
    neg();
    add(sqrtTmp);

    // calculate inverfc(x) = -invphi(x/2)/(sqrt(2))
    mul(SQRT1_2);
    if (sign>0)
      neg();
  }

//*****************************************************************************
// Calendar conversions taken from http://www.fourmilab.ch/documents/calendar/

  private static int floorDiv(int a, int b) {
    if (a>=0)
      return a/b;
    return -((-a+b-1)/b);
  }
  private static int floorMod(int a, int b) {
    if (a>=0)
      return a%b;
    return a+((-a+b-1)/b)*b;
  }

  private static boolean leap_gregorian(int year) {
    return ((year % 4) == 0) &&
      (!(((year % 100) == 0) && ((year % 400) != 0)));
  }

  //  GREGORIAN_TO_JD -- Determine Julian day number from Gregorian
  //  calendar date -- Except that we use 1/1-0 as day 0
  private static int gregorian_to_jd(int year, int month, int day) {
    return ((366 - 1) +
            (365 * (year - 1)) +
            (floorDiv(year - 1, 4)) +
            (-floorDiv(year - 1, 100)) +
            (floorDiv(year - 1, 400)) +
            ((((367 * month) - 362) / 12) +
             ((month <= 2) ? 0 : (leap_gregorian(year) ? -1 : -2)) + day));
  }

  //  JD_TO_GREGORIAN -- Calculate Gregorian calendar date from Julian
  //  day -- Except that we use 1/1-0 as day 0
  private static int jd_to_gregorian(int jd) {
    int wjd, depoch, quadricent, dqc, cent, dcent, quad, dquad,
        yindex, year, yearday, leapadj, month, day;

    wjd = jd;
    depoch = wjd - 366;
    quadricent = floorDiv(depoch, 146097);
    dqc = floorMod(depoch, 146097);
    cent = floorDiv(dqc, 36524);
    dcent = floorMod(dqc, 36524);
    quad = floorDiv(dcent, 1461);
    dquad = floorMod(dcent, 1461);
    yindex = floorDiv(dquad, 365);
    year = (quadricent * 400) + (cent * 100) + (quad * 4) + yindex;
    if (!((cent == 4) || (yindex == 4)))
      year++;
    yearday = wjd - gregorian_to_jd(year, 1, 1);
    leapadj = ((wjd < gregorian_to_jd(year, 3, 1)) ? 0
               : (leap_gregorian(year) ? 1 : 2));
    month = floorDiv(((yearday + leapadj) * 12) + 373, 367);
    day = (wjd - gregorian_to_jd(year, month, 1)) + 1;

    return (year*100+month)*100+day;
  }

  /**
   * Converts this <code>Real</code> from "hours" to "days, hours,
   * minutes and seconds".
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>The format converted to is encoded into the digits of the
   * number (in decimal form):
   * "<code>DDDDhh.mmss</code>". Here "<code>DDDD</code>," is number
   * of days, "<code>hh</code>" is hours (0-23), "<code>mm</code>" is
   * minutes (0-59) and "<code>ss</code>" is seconds
   * (0-59). Additional digits represent fractions of a second.
   *
   * <p>If the number of hours of the input is greater or equal to
   * 8784 (number of hours in year <code>0</code>), the format
   * converted to is instead "<code>YYYYMMDDhh.mmss</code>". Here
   * "<code>YYYY</code>" is the number of years since the imaginary
   * year <code>0</code> in the Gregorian calendar, extrapolated back
   * from year 1582. "<code>MM</code>" is the month (1-12) and
   * "<code>DD</code>" is the day of the month (1-31). See a thorough
   * discussion of date calculations <a
   * href="http://midp-calc.sourceforge.net/Calc.html#DateNote">here</a>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     ?
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     19
   * </td></tr></table>
   */
  public void toDHMS() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    boolean negative = (this.sign!=0);
    abs();

    int D,m;
    long h;
    h = toLong();
    frac();
    tmp1.assign(60);
    mul(tmp1);
    m = toInteger();
    frac();
    mul(tmp1);

    // MAGIC ROUNDING: Check if we are 2**-16 second short of a whole minute
    // i.e. "seconds" > 59.999985
    { tmp2.mantissa = ONE.mantissa; tmp2.exponent = ONE.exponent; tmp2.sign = ONE.sign; };
    tmp2.scalbn(-16);
    add(tmp2);
    if (this.compare(tmp1) >= 0) {
      // Yes. So set zero seconds instead and carry over to minutes and hours
      { this.mantissa = ZERO.mantissa; this.exponent = ZERO.exponent; this.sign = ZERO.sign; };
      m++;
      if (m >= 60) {
        m -= 60;
        h++;
      }
      // Phew! That was close. From now on it is integer arithmetic...
    } else {
      // Nope. So try to undo the damage...
      sub(tmp2);
    }

    D = (int)(h/24);
    h %= 24;
    if (D >= 366)
      D = jd_to_gregorian(D);

    add(m*100);
    div(10000);
    tmp1.assign(D*100L+h);
    add(tmp1);

    if (negative)
      neg();
  }

  /**
   * Converts this <code>Real</code> from "days, hours, minutes and
   * seconds" to "hours".
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p>The format converted from is encoded into the digits of the
   * number (in decimal form):
   * "<code>DDDDhh.mmss</code>". Here "<code>DDDD</code>" is number of
   * days, "<code>hh</code>" is hours (0-23), "<code>mm</code>" is
   * minutes (0-59) and "<code>ss</code>" is seconds
   * (0-59). Additional digits represent fractions of a second.
   *
   * <p>If the number of days in the input is greater than or equal to
   * 10000, the format converted from is instead
   * "<code>YYYYMMDDhh.mmss</code>". Here "<code>YYYY</code>" is the
   * number of years since the imaginary year <code>0</code> in the
   * Gregorian calendar, extrapolated back from year
   * 1582. "<code>MM</code>" is the month (1-12) and
   * "<code>DD</code>" is the day of the month (1-31). If month or day
   * is 0 it is treated as 1. See a thorough discussion of date
   * calculations <a
   * href="http://midp-calc.sourceforge.net/Calc.html#DateNote">here</a>.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     ?
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     19
   * </td></tr></table>
   */
  public void fromDHMS() {
    if (!(this.exponent >= 0 && this.mantissa != 0))
      return;
    boolean negative = (this.sign!=0);
    abs();

    int Y,M,D,m;
    long h;
    h = toLong();
    frac();
    tmp1.assign(100);
    mul(tmp1);
    m = toInteger();
    frac();
    mul(tmp1);

    // MAGIC ROUNDING: Check if we are 2**-10 second short of 100 seconds
    // i.e. "seconds" > 99.999
    { tmp2.mantissa = ONE.mantissa; tmp2.exponent = ONE.exponent; tmp2.sign = ONE.sign; };
    tmp2.scalbn(-10);
    add(tmp2);
    if (this.compare(tmp1) >= 0) {
      // Yes. So set zero seconds instead and carry over to minutes and hours
      { this.mantissa = ZERO.mantissa; this.exponent = ZERO.exponent; this.sign = ZERO.sign; };
      m++;
      if (m >= 100) {
        m -= 100;
        h++;
      }
      // Phew! That was close. From now on it is integer arithmetic...
    } else {
      // Nope. So try to undo the damage...
      sub(tmp2);
    }

    D = (int)(h/100);
    h %= 100;
    if (D>=10000) {
      M = D/100;
      D %= 100;
      if (D==0) D=1;
      Y = M/100;
      M %= 100;
      if (M==0) M=1;
      D = gregorian_to_jd(Y,M,D);
    }

    add(m*60);
    div(3600);
    tmp1.assign(D*24L+h);
    add(tmp1);

    if (negative)
      neg();
  }

  /**
   * Assigns this <code>Real</code> the current time. The time is
   * encoded into the digits of the number (in decimal form), using the
   * format "<code>hh.mmss</code>", where "<code>hh</code>" is hours,
   * "<code>mm</code>" is minutes and "code>ss</code>" is seconds.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     � ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     8.9
   * </td></tr></table>
   */
  public void time() {
    long now = System.currentTimeMillis();
    int h,m,s;
    now /= 1000;
    s = (int)(now % 60);
    now /= 60;
    m = (int)(now % 60);
    now /= 60;
    h = (int)(now % 24);
    assign((h*100+m)*100+s);
    div(10000);
  }

  /**
   * Assigns this <code>Real</code> the current date. The date is
   * encoded into the digits of the number (in decimal form), using
   * the format "<code>YYYYMMDD00</code>", where "<code>YYYY</code>"
   * is the year, "<code>MM</code>" is the month (1-12) and
   * "<code>DD</code>" is the day of the month (1-31). The
   * "<code>00</code>" in this format is a sort of padding to make it
   * compatible with the format used by {@link #toDHMS()} and {@link
   * #fromDHMS()}.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <i>none</i>
   * </td></tr><tr><td><i>Error&nbsp;bound:</i></td><td>
   *     0 ULPs
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     30
   * </td></tr></table>
   */
  public void date() {
    long now = System.currentTimeMillis();
    now /= 86400000; // days
    now *= 24; // hours
    assign(now);
    add(719528*24); // 1970-01-01 era
    toDHMS();
  }

//*****************************************************************************

  /**
   * The seed of the first 64-bit CRC generator of the random
   * routine. Set this value to control the generated sequence of random
   * numbers. Should never be set to 0. See {@link #random()}.
   */
  public static long randSeedA = 0x6487ed5110b4611aL;// Something to start with
  /**
   * The seed of the second 64-bit CRC generator of the random
   * routine. Set this value to control the generated sequence of random
   * numbers. Should never be set to 0. See {@link #random()}.
   */
  public static long randSeedB = 0x56fc2a2c515da54dL;// (mantissas of pi and e)

  // 64 Bit CRC Generators
  //
  // The generators used here are not cryptographically secure, but
  // two weak generators are combined into one strong generator by
  // skipping bits from one generator whenever the other generator
  // produces a 0-bit.
  private static void advanceBit() {
    randSeedA = (randSeedA<<1)^(randSeedA<0?0x1b:0);
    randSeedB = (randSeedB<<1)^(randSeedB<0?0xb000000000000001L:0);
  }

  // Get next bits from the pseudo-random sequence
  private static long nextBits(int bits) {
    long answer = 0;
    while (bits-- > 0) {
      while (randSeedA >= 0)
        advanceBit();
      answer = (answer<<1) + (randSeedB < 0 ? 1 : 0);
      advanceBit();
    }
    return answer;
  }

  /**
   * Accumulate more randomness into the random number generator, to
   * decrease the predictability of the output from {@link
   * #random()}. The input should contain data with some form of
   * inherent randomness e.g. System.currentTimeMillis().
   *
   * @param seed some extra randomness for the random number generator.
   */
  public static void accumulateRandomness(long seed) {
    randSeedA ^= seed & 0x5555555555555555L;
    randSeedB ^= seed & 0xaaaaaaaaaaaaaaaaL;
    nextBits(63);
  }

  /**
   * Calculates a pseudorandom number in the range [0,&nbsp;1).
   * Replaces the contents of this <code>Real</code> with the result.
   *
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this = Math.{@link Math#random() random}();</code>
   * </td></tr><tr><td><i>Approximate&nbsp;error&nbsp;bound:</i></td><td>
   *     -
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     81
   * </td></tr></table>
   */
  public Real random() {
    sign = 0;
    exponent = 0x3fffffff;
    while (nextBits(1) == 0)
      exponent--;
    mantissa = 0x4000000000000000L+nextBits(62);
    return this;
  }

//*****************************************************************************

  private int digit(char a, int base, boolean twosComplement) {
    int digit = -1;
    if (a>='0' && a<='9')
      digit = (int)(a-'0');
    else if (a>='A' && a<='F')
      digit = (int)(a-'A')+10;
    if (digit >= base)
      return -1;
    if (twosComplement)
      digit ^= base-1;
    return digit;
  }

  private void shiftUp(int base) {
    if (base==2)
      scalbn(1);
    else if (base==8)
      scalbn(3);
    else if (base==16)
      scalbn(4);
    else
      mul10();
  }

  private void atof(String a, int base) {
    makeZero();
    int length = a.length();
    int index = 0;
    byte tmpSign = 0;
    boolean compl = false;
    while (index<length && a.charAt(index)==' ')
      index++;
    if (index<length && a.charAt(index)=='-') {
      tmpSign=1;
      index++;
    } else if (index<length && a.charAt(index)=='+') {
      index++;
    } else if (index<length && a.charAt(index)=='/') {
      // Input is twos complemented negative number
      compl=true;
      tmpSign=1;
      index++;
    }
    int d;
    while (index<length && (d=digit(a.charAt(index),base,compl))>=0) {
      shiftUp(base);
      add(d);
      index++;
    }
    if (index<length && (a.charAt(index)=='.' || a.charAt(index)==',')) {
      index++;
      { tmp2.mantissa = ONE.mantissa; tmp2.exponent = ONE.exponent; tmp2.sign = ONE.sign; };
      while (index<length && (d=digit(a.charAt(index),base,compl))>=0) {
        tmp2.shiftUp(base);
        shiftUp(base);
        add(d);
        index++;
      }
      if (compl)
        add(ONE);
      div(tmp2);
    } else if (compl) {
      add(ONE);
    }
    while (index<length && a.charAt(index)==' ')
      index++;
    if (index<length && (a.charAt(index)=='e' || a.charAt(index)=='E')) {
      index++;
      int exp = 0;
      boolean expNeg = false;
      if (index<length && a.charAt(index)=='-') {
        expNeg = true;
        index++;
      } else if (index<length && a.charAt(index)=='+') {
        index++;
      }
      while (index<length && a.charAt(index)>='0' && a.charAt(index)<='9') {
        if (exp < 400000000) // This takes care of overflows and makes inf or 0
          exp = (exp*10)+(int)(a.charAt(index)-'0');
        index++;
      }
      if (expNeg)
        exp = -exp;

      if (base==2)
        scalbn(exp);
      else if (base==8)
        scalbn(exp*3);
      else if (base==16)
        scalbn(exp*4);
      else {
        if (exponent > 300000000 || exponent < -300000000) {
          // Kludge to be able to enter very large and very small numbers
          // without causing over/underflows
          { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
          tmp1.pow(exp/2);
          mul(tmp1);
          exp -= exp/2;
        }
        { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
        tmp1.pow(exp);
        mul(tmp1);
      }
    }
    sign = tmpSign;
  }

//*****************************************************************************

  private void normalizeBCD() {
    if (mantissa == 0) {
      exponent = 0;
      return;
    }
    int carry = 0;
    for (int i=0; i<64; i+=4) {
      int d = (int)((mantissa>>>i)&0xf) + carry;
      carry = 0;
      if (d >= 10) {
        d -= 10;
        carry = 1;
      }
      mantissa &= ~(0xfL<<i);
      mantissa += (long)d<<i;
    }
    if (carry != 0) {
      if ((int)(mantissa&0xf)>=5)
        mantissa += 0x10; // Rounding, may be inaccurate
      mantissa >>>= 4;
      mantissa += 1L<<60;
      exponent++;
      if ((int)(mantissa&0xf)>=10) {
        // Oh, no, not again!
        normalizeBCD();
      }
    }
    while ((mantissa >>> 60) == 0) {
      mantissa <<= 4;
      exponent--;
    }
  }

  private void toBCD() {
    // Convert normal nonzero finite Real to BCD format, represented by 16
    // 4-bit decimal digits in the mantissa, and exponent as a power of 10
    { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
    tmp1.abs();
    { tmp2.mantissa = tmp1.mantissa; tmp2.exponent = tmp1.exponent; tmp2.sign = tmp1.sign; };
    exponent = tmp1.lowPow10();
    if (exponent > 300000000 || exponent < -300000000) {
      // Kludge to be able to print very large and very small numbers
      // without causing over/underflows
      { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
      tmp1.pow(exponent/2);
      tmp2.div(tmp1); // So, divide twice by not-so-extreme numbers
      { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
      tmp1.pow(exponent-(exponent/2));
    }
    tmp2.div(tmp1);
    mantissa = 0;
    for (int i=60; i>=0; i-=4) {
      { tmp1.mantissa = tmp2.mantissa; tmp1.exponent = tmp2.exponent; tmp1.sign = tmp2.sign; };
      tmp1.floor();
      mantissa += (long)tmp1.toInteger()<<i;
      tmp2.sub(tmp1);
      tmp2.mul10();
    }
    if (tmp2.compare(FIVE) >= 0)
      mantissa++; // Rounding
    normalizeBCD();
  }

  private boolean testCarryWhenRounded(int base,int bitsPerDigit,int precision)
  {
    if (mantissa >= 0)
      return false; // msb not set -> carry cannot occur
    int shift = 64-precision*bitsPerDigit;
    if (shift<=0)
      return false; // too high precision -> no rounding necessary
    { tmp5.mantissa = this.mantissa; tmp5.exponent = this.exponent; tmp5.sign = this.sign; };
    if (base==10) {
      tmp5.mantissa += (5L<<(shift-4));
      tmp5.normalizeBCD();
    } else {
      tmp5.mantissa += (1L<<(shift-1));
    }
    if (tmp5.mantissa >= 0) {
      // msb has changed from 1 to 0 -> we have carry
      // now modify the original number
      mantissa = (1L<<(64-bitsPerDigit));
      exponent++;
      return true;
    }
    return false;
  }

  private void round(int base, int bitsPerDigit, int precision) {
    if ((this.exponent == 0 && this.mantissa == 0))
      return;
    int shift = 64-precision*bitsPerDigit;
    if (shift<=0)
      return; // too high precision -> no rounding necessary
    long orgMantissa = mantissa;
    if (base==10) {
      mantissa += (5L<<(shift-4));
      normalizeBCD();
    } else {
      mantissa += (1L<<(shift-1));
    }
    if (orgMantissa < 0 && mantissa >= 0) {
      // msb has changed from 1 to 0 -> we have carry
      mantissa = (1L<<(64-bitsPerDigit));
      exponent++;
    }
  }

  /**
   * The number format used to convert <code>Real</code> values to
   * <code>String</code> using {@link Real#toString(Real.NumberFormat)
   * Real.toString()}. The default number format uses base-10, maximum
   * precision, removal of trailing zeros and '.' as radix point.
   *
   * <p>Note that the fields of the <code>NumberFormat</code> is not
   * protected in any way, the user is responsible for setting the
   * correct values to get a correct result.
   */
  public static class NumberFormat
  {
    /**
     * The number base of the conversion. The default value is 10,
     * valid options are 2, 8, 10 and 16. See {@link Real#and(Real)
     * Real.and()} for an explanation of the interpretation of a
     * <code>Real</code> in base 2, 8 and 16.
     *
     * <p>Negative numbers output in base-2, base-8 and base-16 are
     * shown in two's complement form. This form guarantees that a
     * negative number starts with at least one digit that is the
     * maximum digit for that base, i.e. '1', '7', and 'F',
     * respectively. A positive number is guaranteed to start with at
     * least one '0'. Both positive and negative numbers are extended
     * to the left using this digit, until {@link #maxwidth} is
     * reached.
     */
    public int base = 10;

    /**
     * Maximum width of the converted string. The default value is 30.
     * If the conversion of a <code>Real</code> with a given {@link
     * #precision} would produce a string wider than
     * <code>maxwidth</code>, <code>precision</code> is reduced until
     * the number fits within the given width. If
     * <code>maxwidth</code> is too small to hold the number with its
     * sign, exponent and a <code>precision</code> of 1 digit, the
     * string may become wider than <code>maxwidth</code>.
     *
     * <p>If <code>align</code> is set to anything but
     * <code>ALIGN_NONE</code> and the converted string is shorter
     * than <code>maxwidth</code>, the resulting string is padded with
     * spaces to the specified width according to the alignment.
     */
    public int maxwidth = 30;

    /**
     * The precision, or number of digits after the radix point in the
     * converted string when using the <i>FIX</i>, <i>SCI</i> or
     * <i>ENG</i> format (see {@link #fse}). The default value is 16,
     * valid values are 0-16 for base-10 and base-16 conversion, 0-21
     * for base-8 conversion, and 0-63 for base-2 conversion.
     *
     * <p>The <code>precision</code> may be reduced to make the number
     * fit within {@link #maxwidth}. The <code>precision</code> is
     * also reduced if it is set higher than the actual numbers of
     * significant digits in a <code>Real</code>. When
     * <code>fse</code> is set to <code>FSE_NONE</code>, i.e. "normal"
     * output, the precision is always at maximum, but trailing zeros
     * are removed.
     */
    public int precision = 16;

    /**
     * The special output formats <i>FIX</i>, <i>SCI</i> or <i>ENG</i>
     * are enabled with this field. The default value is
     * <code>FSE_NONE</code>. Valid options are listed below.
     *
     * <p>Numbers are output in one of two main forms, according to
     * this setting. The normal form has an optional sign, one or more
     * digits before the radix point, and zero or more digits after the
     * radix point, for example like this:<br>
     * <code>&nbsp;&nbsp;&nbsp;3.14159</code><br>
     * The exponent form is like the normal form, followed by an
     * exponent marker 'e', an optional sign and one or more exponent
     * digits, for example like this:<br>
     * <code>&nbsp;&nbsp;&nbsp;-3.4753e-13</code>
     *
     * <p><dl>
     *   <dt>{@link #FSE_NONE}
     *   <dd>Normal output. Numbers are output with maximum precision,
     *   trailing zeros are removed. The format is changed to
     *   exponent form if the number is larger than the number of
     *   significant digits allows, or if the resulting string would
     *   exceed <code>maxwidth</code> without the exponent form.
     *
     *   <dt>{@link #FSE_FIX}
     *   <dd>Like normal output, but the numbers are output with a
     *   fixed number of digits after the radix point, according to
     *   {@link #precision}. Trailing zeros are not removed.
     *
     *   <dt>{@link #FSE_SCI}
     *   <dd>The numbers are always output in the exponent form, with
     *   one digit before the radix point, and a fixed number of
     *   digits after the radix point, according to
     *   <code>precision</code>. Trailing zeros are not removed.
     *
     *   <dt>{@link #FSE_ENG}
     *   <dd>Like the <i>SCI</i> format, but the output shows one to
     *   three digits before the radix point, so that the exponent is
     *   always divisible by 3.
     * </dl>
     */
    public int fse = FSE_NONE;

    /**
     * The character used as the radix point. The default value is
     * <code>'.'</code>. Theoretcally any character that does not
     * otherwise occur in the output can be used, such as
     * <code>','</code>.
     *
     * <p>Note that setting this to anything but <code>'.'</code> and
     * <code>','</code> is not supported by any conversion method from
     * <code>String</code> back to <code>Real</code>.
     */
    public char point = '.';

    /**
     * Set to <code>true</code> to remove the radix point if this is
     * the last character in the converted string. This is the
     * default.
     */
    public boolean removePoint = true;

    /**
     * The character used as the thousands separator. The default
     * value is the character code <code>0</code>, which disables
     * thousands-separation. Theoretcally any character that does not
     * otherwise occur in the output can be used, such as
     * <code>','</code> or <code>' '</code>.
     *
     * <p>When <code>thousand!=0</code>, this character is inserted
     * between every 3rd digit to the left of the radix point in
     * base-10 conversion. In base-16 conversion, the separator is
     * inserted between every 4th digit, and in base-2 conversion the
     * separator is inserted between every 8th digit. In base-8
     * conversion, no separator is ever inserted.
     *
     * <p>Note that tousands separators are not supported by any
     * conversion method from <code>String</code> back to
     * <code>Real</code>, so use of a thousands separator is meant
     * only for the presentation of numbers.
     */
    public char thousand = 0;

    /**
     * The alignment of the output string within a field of {@link
     * #maxwidth} characters. The default value is
     * <code>ALIGN_NONE</code>. Valid options are defined as follows:
     *
     * <p><dl>
     *   <dt>{@link #ALIGN_NONE}
     *   <dd>The resulting string is not padded with spaces.
     *
     *   <dt>{@link #ALIGN_LEFT}
     *   <dd>The resulting string is padded with spaces on the right side
     *   until a width of <code>maxwidth</code> is reached, making the
     *   number left-aligned within the field.
     *
     *   <dt>{@link #ALIGN_RIGHT}
     *   <dd>The resulting string is padded with spaces on the left side
     *   until a width of <code>maxwidth</code> is reached, making the
     *   number right-aligned within the field.
     *
     *   <dt>{@link #ALIGN_CENTER}
     *   <dd>The resulting string is padded with spaces on both sides
     *   until a width of <code>maxwidth</code> is reached, making the
     *   number center-aligned within the field.
     * </dl>
     */
    public int align = ALIGN_NONE;

    /** Normal output {@linkplain #fse format} */
    public static final int FSE_NONE = 0;
    /** <i>FIX</i> output {@linkplain #fse format} */
    public static final int FSE_FIX = 1;
    /** <i>SCI</i> output {@linkplain #fse format} */
    public static final int FSE_SCI = 2;
    /** <i>ENG</i> output {@linkplain #fse format} */
    public static final int FSE_ENG = 3;

    /** No {@linkplain #align alignment} */
    public static final int ALIGN_NONE = 0;
    /** Left {@linkplain #align alignment} */
    public static final int ALIGN_LEFT = 1;
    /** Right {@linkplain #align alignment} */
    public static final int ALIGN_RIGHT = 2;
    /** Center {@linkplain #align alignment} */
    public static final int ALIGN_CENTER = 3;
  }

  private String align(StringBuffer s, NumberFormat format) {
    if (format.align == NumberFormat.ALIGN_LEFT) {
      while (s.length()<format.maxwidth)
        s.append(' ');
    } else if (format.align == NumberFormat.ALIGN_RIGHT) {
      while (s.length()<format.maxwidth)
        s.insert(0,' ');
    } else if (format.align == NumberFormat.ALIGN_CENTER) {
      while (s.length()<format.maxwidth) {
        s.append(' ');
        if (s.length()<format.maxwidth)
          s.insert(0,' ');
      }
    }
    return s.toString();
  }

  private static StringBuffer ftoaBuf = new StringBuffer(40);
  private static StringBuffer ftoaExp = new StringBuffer(15);
  /**
   * This string holds the only valid characters to use in hexadecimal
   * numbers. Equals <code>"0123456789ABCDEF"</code>.
   * See {@link #assign(String,int)}.
   */
  public static final String hexChar = "0123456789ABCDEF";

  private String ftoa(NumberFormat format) {
    ftoaBuf.setLength(0);

    if ((this.exponent < 0 && this.mantissa != 0)) {
      ftoaBuf.append("nan");
      return align(ftoaBuf,format);
    }
    if ((this.exponent < 0 && this.mantissa == 0)) {
      ftoaBuf.append((this.sign!=0) ? "-inf":"inf");
      return align(ftoaBuf,format);
    }

    int bitsPerDigit,digitsPerThousand;
    switch (format.base) {
      case 2:
        bitsPerDigit = 1;
        digitsPerThousand = 8;
        break;
      case 8:
        bitsPerDigit = 3;
        digitsPerThousand = 1000; // Disable thousands separator
        break;
      case 16:
        bitsPerDigit = 4;
        digitsPerThousand = 4;
        break;
      case 10:
      default:
        bitsPerDigit = 4;
        digitsPerThousand = 3;
        break;
    }
    if (format.thousand == 0)
      digitsPerThousand = 1000; // Disable thousands separator

    int accurateBits = 64;
    { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; };
    if ((this.exponent == 0 && this.mantissa == 0)) {
      tmp4.exponent = 0;
      if (format.base != 10)
        tmp4.sign = 0; // Two's complement cannot display -0
    } else {
      if (format.base == 10) {
        tmp4.toBCD();
      } else {
        if ((tmp4.sign!=0)) {
          tmp4.mantissa = -tmp4.mantissa;
          if (((tmp4.mantissa >> 62)&3) == 3) {
            tmp4.mantissa <<= 1;
            tmp4.exponent--;
            accurateBits--; // ?
          }
        }
        tmp4.exponent -= 0x40000000-1;
        int shift = bitsPerDigit-1 - floorMod(tmp4.exponent,bitsPerDigit);
        tmp4.exponent = floorDiv(tmp4.exponent,bitsPerDigit);
        if (shift == bitsPerDigit-1) {
          // More accurate to shift up instead
          tmp4.mantissa <<= 1;
          tmp4.exponent--;
          accurateBits--;
        }
        else if (shift>0) {
          tmp4.mantissa = (tmp4.mantissa+(1L<<(shift-1)))>>>shift;
          if ((tmp4.sign!=0)) // Need to fill in some 1's at the top
            tmp4.mantissa |= 0x8000000000000000L>>(shift-1); // >> not >>>
        }
      }
    }
    int accurateDigits = (accurateBits+bitsPerDigit-1)/bitsPerDigit;

    int precision;
    int pointPos = 0;
    do
    {
      int width = format.maxwidth-1; // subtract 1 for decimal point
      int prefix = 0;
      if (format.base != 10)
        prefix = 1; // want room for at least one "0" or "f/7/1"
      else if ((tmp4.sign!=0))
        width--; // subtract 1 for sign

      boolean useExp = false;
      switch (format.fse) {
        case NumberFormat.FSE_SCI:
          precision = format.precision+1;
          useExp = true;
          break;
        case NumberFormat.FSE_ENG:
          pointPos = floorMod(tmp4.exponent,3);
          precision = format.precision+1+pointPos;
          useExp = true;
          break;
        case NumberFormat.FSE_FIX:
        case NumberFormat.FSE_NONE:
        default:
          precision = 1000;
          if (format.fse == NumberFormat.FSE_FIX)
            precision = format.precision+1;
          if (tmp4.exponent+1 >
                width-(tmp4.exponent+prefix)/digitsPerThousand-prefix+
                (format.removePoint ? 1:0) ||
              tmp4.exponent+1 > accurateDigits ||
              -tmp4.exponent >= width ||
              -tmp4.exponent >= precision)
          {
            useExp = true;
          } else {
            pointPos = tmp4.exponent;
            precision += tmp4.exponent;
            if (tmp4.exponent > 0)
              width -= (tmp4.exponent+prefix)/digitsPerThousand;
            if (format.removePoint && tmp4.exponent == width-prefix)
              width++; // Add  1 for the decimal point that will be removed
          }
          break;
      }

      if (prefix!=0 && pointPos>=0)
        width -= prefix;

      ftoaExp.setLength(0);
      if (useExp) {
        ftoaExp.append('e');
        ftoaExp.append(tmp4.exponent-pointPos);
        width -= ftoaExp.length();
      }
      if (precision > accurateDigits)
        precision = accurateDigits;
      if (precision > width)
        precision = width;
      if (precision > width+pointPos) // In case of negative pointPos
        precision = width+pointPos;
      if (precision <= 0)
        precision = 1;
    }
    while (tmp4.testCarryWhenRounded(format.base,bitsPerDigit,precision));

    tmp4.round(format.base,bitsPerDigit,precision);

    // Start generating the string. First the sign
    if ((tmp4.sign!=0) && format.base == 10)
      ftoaBuf.append('-');

    // Save pointPos for hex/oct/bin prefixing with thousands-sep
    int pointPos2 = pointPos;

    // Add leading zeros (or f/7/1)
    char prefixChar = (format.base==10 || (tmp4.sign==0)) ? '0' :
                      hexChar.charAt(format.base-1);
    if (pointPos < 0) {
      ftoaBuf.append(prefixChar);
      ftoaBuf.append(format.point);
      while (pointPos < -1) {
        ftoaBuf.append(prefixChar);
        pointPos++;
      }
    }

    // Add fractional part
    while (precision > 0) {
      ftoaBuf.append(hexChar.charAt((int)(tmp4.mantissa>>>(64-bitsPerDigit))));
      tmp4.mantissa <<= bitsPerDigit;
      if (pointPos>0 && pointPos%digitsPerThousand==0)
        ftoaBuf.append(format.thousand);
      if (pointPos == 0)
        ftoaBuf.append(format.point);
      precision--;
      pointPos--;
    }

    if (format.fse == NumberFormat.FSE_NONE) {
      // Remove trailing zeros
      while (ftoaBuf.charAt(ftoaBuf.length()-1) == '0')
        ftoaBuf.setLength(ftoaBuf.length()-1);
    }
    if (format.removePoint) {
      // Remove trailing point
      if (ftoaBuf.charAt(ftoaBuf.length()-1) == format.point)
        ftoaBuf.setLength(ftoaBuf.length()-1);
    }

    // Add exponent
    ftoaBuf.append(ftoaExp);

    // In case hex/oct/bin number, prefix with 0's or f/7/1's
    if (format.base!=10) {
      while (ftoaBuf.length()<format.maxwidth) {
        pointPos2++;
        if (pointPos2>0 && pointPos2%digitsPerThousand==0)
          ftoaBuf.insert(0,format.thousand);
        if (ftoaBuf.length()<format.maxwidth)
          ftoaBuf.insert(0,prefixChar);
      }
      if (ftoaBuf.charAt(0) == format.thousand)
        ftoaBuf.deleteCharAt(0);
    }

    return align(ftoaBuf,format);
  }

  private static NumberFormat tmpFormat = new NumberFormat();

  /**
   * Converts this <code>Real</code> to a <code>String</code> using
   * the default <code>NumberFormat</code>.
   *
   * <p>See {@link Real.NumberFormat NumberFormat} for a description
   * of the default way that numbers are formatted.
   * 
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td><td>
   *     <code>this.toString()
   * </td></tr><tr><td><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;
   * </i></td><td>
   *     150
   * </td></tr></table>
   *
   * @return a <code>String</code> representation of this <code>Real</code>.
   */
  public String toString() {
    tmpFormat.base = 10;
    return ftoa(tmpFormat);
  }

  /**
   * Converts this <code>Real</code> to a <code>String</code> using
   * the default <code>NumberFormat</code> with <code>base</code> set
   * according to the argument.
   *
   * <p>See {@link Real.NumberFormat NumberFormat} for a description
   * of the default way that numbers are formatted.
   * 
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td>
   * <td colspan="2">
   *     <code>this.toString()  // Works only for base-10</code>
   * </td></tr><tr><td rowspan="4" valign="top"><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;</i>
   * </td><td width="1%">base-2</td><td>
   *     100
   * </td></tr><tr><td>base-8</td><td>
   *     110
   * </td></tr><tr><td>base-10</td><td>
   *     150
   * </td></tr><tr><td>base-16&nbsp;&nbsp;</td><td>
   *     120
   * </td></tr></table>
   *
   * @param base the base for the conversion. Valid base values are
   *     2, 8, 10 and 16.
   * @return a <code>String</code> representation of this <code>Real</code>.
   */
  public String toString(int base) {
    tmpFormat.base = base;
    return ftoa(tmpFormat);
  }

  /**
   * Converts this <code>Real</code> to a <code>String</code> using
   * the given <code>NumberFormat</code>.
   *
   * <p>See {@link Real.NumberFormat NumberFormat} for a description of the
   * various ways that numbers may be formatted.
   * 
   * <p><table border="1" width="100%" cellpadding="3" cellspacing="0"
   * bgcolor="#e8d0ff"><tr><td width="1%"><i>
   * Equivalent&nbsp;</i><code>double</code><i>&nbsp;code:</i></td>
   * <td colspan="2">
   *     <code>String.format("%...g",this);  // Works only for base-10</code>
   * </td></tr><tr><td rowspan="4" valign="top"><i>
   * Execution&nbsp;time&nbsp;relative&nbsp;to&nbsp;add:&nbsp;&nbsp;</i>
   * </td><td width="1%">base-2</td><td>
   *     100
   * </td></tr><tr><td>base-8</td><td>
   *     110
   * </td></tr><tr><td>base-10</td><td>
   *     150
   * </td></tr><tr><td>base-16&nbsp;&nbsp;</td><td>
   *     120
   * </td></tr></table>
   *
   * @param format the number format to use in the conversion.
   * @return a <code>String</code> representation of this <code>Real</code>.
   */
  public String toString(NumberFormat format) {
    return ftoa(format);
  }
}
